两点边值问题求解
原问题求解思路一阶收敛速度下的迭代思想方法一阶实现代码二阶收敛速度下的迭代思想方法二阶实现代码
实验结果迭代结果误差收敛速度
小结
原问题
−
u
′
′
(
x
)
+
u
(
x
)
=
e
x
(
sin
x
−
2
cos
x
)
,
0
⩽
x
⩽
π
-u^{\prime \prime}(x)+u(x)=e^{x}(\sin x-2 \cos x), \quad 0 \leqslant x \leqslant \pi
−u′′(x)+u(x)=ex(sinx−2cosx),0⩽x⩽π
u
(
0
)
=
0
,
u
(
π
)
=
0
u(0)=0, \quad u(\pi)=0
u(0)=0,u(π)=0
求解思路
一阶收敛速度下的迭代思想方法
(
1
0
−
1
2
+
h
2
q
(
x
1
)
−
1
⋱
⋱
⋱
−
1
2
+
h
2
q
(
x
m
−
1
)
−
1
0
1
)
(
u
0
u
1
⋮
u
m
−
1
u
m
)
\left(\begin{array}{ccccc}1 & 0 & & & \\ -1 & 2+h^{2} q\left(x_{1}\right) & -1 & & \\ & \ddots & \ddots & \ddots & \\ & & -1 & 2+h^{2} q\left(x_{m-1}\right) & -1 \\ & & & 0 & 1\end{array}\right)\left(\begin{array}{c}u_{0} \\ u_{1} \\ \vdots \\ u_{m-1} \\ u_{m}\end{array}\right)
⎝⎜⎜⎜⎜⎛1−102+h2q(x1)⋱−1⋱−1⋱2+h2q(xm−1)0−11⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎛u0u1⋮um−1um⎠⎟⎟⎟⎟⎟⎞
=
(
α
h
2
f
(
x
1
)
⋮
h
2
f
(
x
m
−
1
)
β
)
=\left(\begin{array}{c}\alpha \\ h^{2} f\left(x_{1}\right) \\ \vdots \\ h^{2} f\left(x_{m-1}\right) \\ \beta\end{array}\right)
=⎝⎜⎜⎜⎜⎜⎛αh2f(x1)⋮h2f(xm−1)β⎠⎟⎟⎟⎟⎟⎞ 综合上述问题易得，h=π/n、q(x)=1。
一阶实现代码
%n为步数，h为步长 Au
=E，u为数值解，u1解析解。
%power乘冥 diag短矩阵，用法diag（V，k）将向量V放入到第k条对角线上。
%步长取pi
/10
n
=10;
h
=pi
/n
;
%构造系数矩阵A三对角矩阵，由a、b、c三条对角角线组成
a
=ones
(n
+1,1);
b
=zeros
(n
,1);
c
=zeros
(n
,1);
for i
=2:n
a
(i
,1)=2+power
(h
,2);
end
for i
=1:n
-1
b
(i
,1)=-1;
end
for i
=2:n
c
(i
,1)=-1;
end
B
=diag
(a
);
C
=diag
(b
,-1);
D
=diag
(c
,1);
A
=B
+C
+D
;
%构造系数向量组E
%迭代函数f（x）
=exp
(x
).*(sin
(x
)-2cos
(x
))
d
=zeros
(n
+1,1);
for i
=2:n
d
(i
,1)=power
(h
,2)*exp
(h
*(i
-1)).*(sin
(h
*(i
-1))-2*cos
(h
*(i
-1)));
end
%求数值解向量u
u1
=A\d
;
%步长取pi
/20
n
=20;
h
=pi
/n
;
%构造系数矩阵A三对角矩阵，由a、b、c三条对角角线组成
a
=ones
(n
+1,1);
b
=zeros
(n
,1);
c
=zeros
(n
,1);
for i
=2:n
a
(i
,1)=2+power
(h
,2);
end
for i
=1:n
-1
b
(i
,1)=-1;
end
for i
=2:n
c
(i
,1)=-1;
end
B
=diag
(a
);
C
=diag
(b
,-1);
D
=diag
(c
,1);
A
=B
+C
+D
;
%构造系数向量组E
%迭代函数f（x）
=exp
(x
).*(sin
(x
)-2cos
(x
))
d
=zeros
(n
+1,1);
for i
=2:n
d
(i
,1)=power
(h
,2)*exp
(h
*(i
-1)).*(sin
(h
*(i
-1))-2*cos
(h
*(i
-1)));
end
%求数值解向量u
u2
=A\d
;
%步长取pi
/40
n
=40;
h
=pi
/n
;
%构造系数矩阵A三对角矩阵，由a、b、c三条对角角线组成
a
=ones
(n
+1,1);
b
=zeros
(n
,1);
c
=zeros
(n
,1);
for i
=2:n
a
(i
,1)=2+power
(h
,2);
end
for i
=1:n
-1
b
(i
,1)=-1;
end
for i
=2:n
c
(i
,1)=-1;
end
B
=diag
(a
);
C
=diag
(b
,-1);
D
=diag
(c
,1);
A
=B
+C
+D
;
%构造系数向量组E
%迭代函数f（x）
=exp
(x
).*(sin
(x
)-2cos
(x
))
d
=zeros
(n
+1,1);
for i
=2:n
d
(i
,1)=power
(h
,2)*exp
(h
*(i
-1)).*(sin
(h
*(i
-1))-2*cos
(h
*(i
-1)));
end
%求数值解向量u
u3
=A\d
;
%步长取pi
/10
n
=80;
h
=pi
/n
;
%构造系数矩阵A三对角矩阵，由a、b、c三条对角角线组成
a
=ones
(n
+1,1);
b
=zeros
(n
,1);
c
=zeros
(n
,1);
for i
=2:n
a
(i
,1)=2+power
(h
,2);
end
for i
=1:n
-1
b
(i
,1)=-1;
end
for i
=2:n
c
(i
,1)=-1;
end
B
=diag
(a
);
C
=diag
(b
,-1);
D
=diag
(c
,1);
A
=B
+C
+D
;
%构造系数向量组E
%迭代函数f（x）
=exp
(x
).*(sin
(x
)-2cos
(x
))
d
=zeros
(n
+1,1);
for i
=2:n
d
(i
,1)=power
(h
,2)*exp
(h
*(i
-1)).*(sin
(h
*(i
-1))-2*cos
(h
*(i
-1)));
end
%求数值解向量u
u4
=A\d
;
%求解析解u
(x
)=exp
(x
)*sin
(x
)
x1
=0:pi
/10:pi
;
x2
=0:pi
/20:pi
;
x3
=0:pi
/40:pi
;
x4
=0:pi
/80:pi
t
=0:0.001:pi
;
y
=exp
(t
).*sin
(t
);
% plot
(x1
,u1
,'*',t
,y
,'r')
% legend
('pi/10','解析解')
plot
(x4
,u4
,'*',t
,y
,'r')
legend
('pi/80','解析解')
二阶收敛速度下的迭代思想方法
(
2
+
h
2
g
(
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
)
\left(\begin{array}{ccccc}2+h^{2} g\left(x_{1}\right) & -1 & & & \\ -1 & 2+h^{2} q\left(x_{2}\right) & -1 & & \\ \ddots & \ddots & \ddots & \\ & & -1 & 2+h^{2} q\left(x_{m-2}\right) & -1 \\ & & & -1 & 2+h^{2} q\left(x_{m-1}\right)\end{array}\right)\left(\begin{array}{c}u_{1} \\ u_{2} \\ \vdots \\ u_{m-2} \\ u_{m-1}\end{array}\right)
⎝⎜⎜⎜⎜⎛2+h2g(x1)−1⋱−12+h2q(x2)⋱−1⋱−12+h2q(xm−2)−1−12+h2q(xm−1)⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎛u1u2⋮um−2um−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{array}{c}h^{2} f\left(x_{1}\right)+\alpha \\ h^{2} f\left(x_{2}\right) \\ \vdots \\ h^{2} f\left(x_{m-2}\right) \\ h^{2} f\left(x_{m-1}\right)+\beta\end{array}\right)
=⎝⎜⎜⎜⎜⎜⎛h2f(x1)+αh2f(x2)⋮h2f(xm−2)h2f(xm−1)+β⎠⎟⎟⎟⎟⎟⎞
综合上述问题易得，h=π/n、q(x)=1。
二阶实现代码
%n为步数，h为步长 Au
=E，u为数值解，u1解析解。
%power乘冥 diag短矩阵，用法diag（V，k）将向量V放入到第k条对角线上。
%步长取pi
/10
n
=80;
h
=pi
/n
;
%构造系数矩阵A三对角矩阵，由a、b、c三条对角角线组成
a
=ones
(n
-1,1);
b
=zeros
(n
-2,1);
c
=zeros
(n
-2,1);
for i
=1:n
-1
a
(i
,1)=2+power
(h
,2);
end
for i
=1:n
-2
b
(i
,1)=-1;
end
for i
=1:n
-2
c
(i
,1)=-1;
end
B
=diag
(a
);
C
=diag
(b
,-1);
D
=diag
(c
,1);
A
=B
+C
+D
;
%构造系数向量组E
%迭代函数f（x）
=exp
(x
).*(sin
(x
)-2cos
(x
))
%apha beta都是
0，不做特别处理
d
=zeros
(n
-1,1);
for i
=1:n
-1
d
(i
,1)=power
(h
,2)*exp
(h
*(i
)).*(sin
(h
*(i
))-2*cos
(h
*(i
)));
end
%求数值解向量u
u1
=A\d
;
u11
=zeros
(n
+1,1);
for i
=2:n
u11
(i
,1)=u1
(i
-1,1);
end
x1
=0:pi
/n
:pi
;
t
=0:0.001:pi
;
y
=exp
(t
).*sin
(t
);
plot
(x1
,u11
,'*',t
,y
,'r')
legend
('pi/80','解析解')
实验结果
迭代结果
“一阶无穷小”差分形式不同步长下的数值解与解析解如下图所示： “二阶无穷小”差分形式不同步长下数值解与解析解如下图所示：
误差收敛速度
以“一阶无穷小”差分为例，（同理得到二阶限于时间精力在此省去）如下表1所示： 不难得到收敛速度为0（
h
2
h^2
h2）
小结
本次实验通过分别用“一二阶无穷小”差分法格式对两点边值问题进行求解，通过本次实验加深了笔者数值分析方面的理论与实践知识，但在进一步的误差分析上，实验仍有待提高。
