简单的差分格式解一个一元二阶常微分方程
对于一个如下的一元二阶常微分方程
{
−
u
′
′
+
q
(
x
)
u
=
f
(
x
)
a
<
x
<
b
u
(
a
)
=
α
u
(
b
)
=
β
\begin{cases} -u''+q\left( x \right) u=f\left( x \right) \,\, a<x<b\\ u\left( a \right) =\alpha \,\,u\left( b \right) =\beta\\ \end{cases}
{−u′′+q(x)u=f(x)a<x<bu(a)=αu(b)=β
其中二阶导数可以换成差分格式
u
′
′
(
x
i
)
=
1
h
2
(
u
(
x
i
−
1
)
−
2
u
(
x
i
)
+
u
(
x
i
+
1
)
)
u''\left( x_i \right) =\frac{1}{h^2}\left( u\left( x_{i-1} \right) -2u\left( x_i \right) +u\left( x_{i+1} \right) \right)
u′′(xi)=h21(u(xi−1)−2u(xi)+u(xi+1))
将
u
(
x
i
)
u(x_i)
u(xi)简写为
u
i
u_i
ui
−
u
i
−
1
h
2
+
(
q
(
x
i
)
+
2
h
2
)
u
i
−
u
i
+
1
h
2
=
f
(
x
i
)
-\frac{u_{i-1}}{h^2}+\left( q\left( x_i \right) +\frac{2}{h^2} \right) u_i-\frac{u_{i+1}}{h^2}=f\left( x_i \right)
−h2ui−1+(q(xi)+h22)ui−h2ui+1=f(xi)
将a到b的区域平均分成m份，解以上的式子可以看做解一个三对角方程组
(
2
+
h
2
q
(
x
1
)
−
1
−
1
2
+
h
2
q
(
x
2
)
−
1
⋱
⋱
⋱
−
1
2
+
h
2
q
(
x
m
−
2
)
−
1
−
1
2
+
h
2
q
(
x
m
−
1
)
)
(
u
1
u
2
⋮
u
m
−
2
u
m
−
1
)
=
(
h
2
f
(
x
1
)
+
α
h
2
f
(
x
2
)
⋮
h
2
f
(
x
m
−
2
)
h
2
f
(
x
m
−
1
)
+
β
)
\left( \begin{matrix} \begin{matrix} 2+h^2q(x_1)& -1& \\ -1& 2+h^2q(x_2)& -1\\ & \ddots& \ddots\\ \end{matrix}& \begin{matrix} & & \\ & & \\ \ddots& & \\ \end{matrix}\\ \begin{matrix} & & \\ & & \\ & & \\ \end{matrix}& \begin{matrix} -1& 2+h^2q(x_{m-2})& -1\\ & -1& 2+h^2q(x_{m-1})\\ & & \\ \end{matrix}\\ \end{matrix} \right) \left( \begin{array}{c} u_1\\ u_2\\ \begin{array}{c} \vdots\\ u_{m-2}\\ \end{array}\\ u_{m-1}\\ \end{array} \right) =\left( \begin{array}{c} h^2f\left( x_1 \right) +\alpha\\ h^2f\left( x_2 \right)\\ \begin{array}{c} \vdots\\ h^2f\left( x_{m-2} \right)\\ \end{array}\\ h^2f\left( x_{m-1} \right) +\beta\\ \end{array} \right)
2+h2q(x1)−1−12+h2q(x2)⋱−1⋱⋱−12+h2q(xm−2)−1−12+h2q(xm−1)
u1u2⋮um−2um−1
=
h2f(x1)+αh2f(x2)⋮h2f(xm−2)h2f(xm−1)+β
这里是左右两边同时乘以
h
2
h^2
h2之后展开成的矩阵形式，取
f
(
x
)
=
e
x
(
sin
(
x
)
−
2
cos
(
x
)
)
f(x)=e^x(\sin(x)-2\cos(x))
f(x)=ex(sin(x)−2cos(x))
q
(
x
)
=
1
q(x)=1
q(x)=1代码如下
#include<stdio.h>
#include<math.h>
#define M_PI 3.14159265358979323846
#define M 20 //m等分
#define N M-1
/*
a、b、c是三对角矩阵的下对角线、主对角线和上对角线的元素，d是右侧向量，x是结果向量
*/
void thomas(double a[], double b[], double c[], double d[], double x[])
{
double y[N],p[N],q[N];
int i;
p[0]=b[0];
for(i=0;i<N-1;i++){
q[i]=c[i]/p[i];
p[i+1]=b[i+1]-a[i]*q[i];
}
y[0]=d[0]/p[0];
for(i=1;i<N;i++){
y[i]=(d[i]-a[i-1]*y[i-1])/p[i];
}
x[N-1]=y[N-1];
for(i=N-2;i>=0;i--){
x[i]=y[i]-q[i]*x[i+1];
}
}
double f(double x){
return exp(x)*(sin(x)-2*cos(x));
}
double q(double x){
return 1;
}
int main(){
double a=0.0;//下限
double b=M_PI;//上限
double h=(b-a)/M;//步长
double alpha=0.0;
double beta=0.0;
double zhu[N],xia[N-1],shang[N-1],you[N],x[N];
int i;
for(i=0;i<N;i++){
zhu[i]=q(a+i*h+h)*h*h+2;
you[i]=f(a+i*h+h)*h*h;
}
you[0]+=alpha;
you[N-1]+=beta;
for(i=0;i<N-1;i++){
xia[i]=-1;
shang[i]=-1;
}
thomas(xia,zhu,shang,you,x);
for(i=0;i<N;i++){
printf("u[%d]=%lf\n",i+1,x[i]);
}
return 0;
}
这里也包含了追赶法求三对角矩阵的函数
原文地址：https://blog.csdn.net/2301_76273740/article/details/140226328
免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！