一阶格式和二阶格式误差对比

Π/1280
1.0972
3.3330
6.2475
Π/2560
1.0995
3.3373
6.2566
精确解
1.1018
3.3416
6.2637
4Π/5 7.0125 7.1348 7.1957 7.2261 7.2412 7.3264
取不同步长在 x / 5,2 / 5,3 / 5,4 / 5 处解析解和数值解的误差矩
对角矩阵，通过循环对于这个矩阵进行赋值，将 u1, u2 um1 放入
向量 y 中，这样方程右边为 b h2 f ( y)
赋值程序为
clc,clear h=pi/160; x=zeros(1,pi/h+1); for i=1:pi/h+1
x(i)=(i-1)*h; end
y=x(1:pi/h+1);%x1 到 x(m+1)组成的矩阵
即可以化成方程 Au=b，若 A 非奇异则 u=inv(A)*b，则可以得到解。
所得到 u0, u1, u2 um 即为以 h 为网格步长得到的网格结点的数值
解。
（三）程序的思路
首先在 matlab 中建立两个函数文件 f.m 和 q.m，这两个函数文件
分别代表偏微分方程中的 f(x)和 q(x)。根据题目
以 h=pi/160 为例得到的图像为：
我们发现 h 的值取的越小时，解析解和数值解所得到的图像越接近， 为了证实这个我们对数值解的误差进行研究。 我们分别取 h /10, / 20, / 40, / 80, /160 进行研究，这样得到数值
解，我们取 x / 5,2 / 5,3 / 5,4 / 5 进行研究并放在一个矩阵中
接下来的数值解进行比较，以此来判断求解方法是否精确。
A 一阶格式
ui1
2ui h2
ui1
qiui
fi
(u1 u0 ) h, um um1 h ,1 i m
将 h2 乘到该方程的两端可得到：
ui1 (u1
[2 h2q(xi u0) ,um
)]ui um1
ui1 h
h ,1
h2 heπ
f ( xm1) h2 / 2 f
(
xm
)
以 h=pi/10 为例得到的图像为
我们发现 h 的值取的越小时，解析解和数值解所得到的图像越接近， 为了证实这个我们对数值解的误差进行研究。
我们分别取 h /10, / 20, / 40, / 80, /160 进行研究，这样得到数值
1.1047
3.3460
6.2828
Π/40
1.1025
3.3427
6.2685
Π/80
1.1020
3.3419
6.2649
Π/160
1.1018
3.3417
6.2640
精确解
1.1018
3.3416
6.2637
4Π/5 7.5018 7.3184 7.2719 7.2603 7.2573 7.3264
Π/640 0.0092 3.3243 0.0324 0.0607 0.1137 2.003
Π/1280 0.0046 0.0086 0.0162 0.0303 0.0568 2.001
Π/2560 0.0023 0.0043 0.0081 0.0152 0.0284 2.000
取不同步长时所取得数值解的误差曲线
所以当步长缩小为原来的 1/2 时，最大误差约缩小为原来的 1/2.
B 二阶格式的解
二阶格式的方程为
1 h2 / 2 1
u0 h h2 / 2 f (x0 )
1 2 h2 1
u1
h2 f ( x1)
1
2 h2 1
1 um1 1 h2 / 2 um
导数边界值问题一阶格式和二阶格式的误差比较
06140211 陈永鑫
（一）问题
应用差分格式计算如下两点边值问题对比一阶格式和二阶格式误差
u''(x) u(x) ex (sin x 2 cos x)
u'
(0)
1,
u'
(
)
e
（二）求解方法描述 首先该方程有解析解 u( x) ex sin x ，有了这个以后可以与我们
取不同步长时所取得数值解的误差曲线为
所以当步长缩小为原来的 1/2 时，最大误差约缩小为原来的 1/4.
（五）结论
一阶格式的误差为 O(h),二阶格式的误差为 O(h²)
fy=zeros(1,pi/h+1);%f(x1)到 f(x(m+1))组成的矩阵 for i=1:pi/h+1
fy(i)=f(y(i)); end
A=zeros(pi/h+1,pi/h+1);%系数矩阵为 A for i=1:pi/h+1
A(i,i)=2+h*h*q(x(i)); end for i=1:pi/h
2 f ( xi i
) m
因为在这个题目中 q(x)=1， 1, eπ ， 1 0,2 0 所以
将这个线性方程写成矩阵形式：
1 1 1 2 h2
1 1 2 h2
1
u0 h
u1 h2 f (x1)
1 um1 h2 f ( xm1)
1 um heπ
f.m 为
function y=f(x); y=exp(x)*(sin(x)-2*cos(x));
q.m 为
function y=q(x); y=1;
这样当题目中的 f(x)和 q(x)需要发生变化时只要改变 f.m 和 q.m 中的
函数即可，不需要去改变主程序。
在 h 设定以后， m ( ) / h 1，所以系数矩阵为(m-1)*(m-1)的三
a=0:0.001:pi; b=exp(a).*sin(a);
plot(a,b)%绘制解析解的图 hold on
plot(X,Y,'-*')%绘制数值解的图 hold off
legend('u(x)数值解','u(x)解析解');
title('h=pi/160 的数值解和解析解的图像')
（四）程序的结果
阵
取不同步长时部分结点处数值解的误差的绝对值和数值解的最大误
差
h\x
Π/5 2Π/5 3Π/5 4Π/5 E(h)
E(2h)/E(h)
Π/160 0.0373 0.0697 0.1304 0.2438 0.4565 *
Π/320 0.0185 0.0347 0.0649 0.1216 0.2277 2.005
A(i,i+1)=-1; end for i=2:pi/h+1
A(i,i-1)=-1; end A(1,1)=1; A(pi/h+1,pi/h+1)=1; b=(h*h*fy)'; b(1)=-h;
b(pi/h+1)=-exp(pi)*h;%方程右端
求解程序为：
u=inv(A)*b;%求得方程的解
画图程序为：
取不同步长在 x / 5,2 / 5,3 / 5,4 / 5 处解析解和数值解的误差矩
阵
取不同步长时部分结点处数值解的误差的绝对值和数值解的最大误
差
h\x
Π/5 2Π/5 3Π/5 4Π/5 E(h)
E(2h)/E(h)
Π/10 0.0116 0.0173 0.0753 0.2455 0.6041 *
解，我们取 x / 5,2 / 5,3 / 5,4 / 5 进行研究并放在一个矩阵中
又解析解 u( x) ex sin x 可以得到精确值为 [1.1018 3.3416 6.2637 7.2564]
即部分结点的取值表为
h\x
Π/5
2Π/5
3Π/5
Π/10
1.1134
3.3590
6.3390
Π/20
Π/20 0.0029 0.0044 0.0191 0.0620 0.1524 3.9639
Π/40 0.0007 0.0011 0.0048 0.0155 0.0382 3.9895
Π/80 0.0002 0.0003 0.0012 0.0039 0.0096 3.979
Π/160 0.0000 0.0001 0.0003 0.0010 0.0024 4.000
又解析解 u( x) ex sin x 可以得到精确值为 [1.1018 3.3416 6.2637 7.2564]
即部分结点Leabharlann Baidu取值表为
h\x
Π/5
2Π/5
3Π/5
Π/160
1.0644
3.2719
6.1334
Π/320
1.0832
3.3069
6.1988
Π/640
1.0925
3.3243
6.2313