常微分方程数值解实验报告.docx
《常微分方程数值解实验报告.docx》由会员分享，可在线阅读，更多相关《常微分方程数值解实验报告.docx（11页珍藏版）》请在冰豆网上搜索。
常微分方程数值解实验报告
常微分方程数值解实验报告
学院：
数学与信息科学
专业：
信息与计算科学
：
思义
学号：
201216524
课程：
常微分方程数值解
实验一：
常微分方程的数值解法
1、分别用Euler法、改进的Euler法（预报校正格式）和S—K法求解初值问题。
（h=0.1）并与真解作比较。
1.1实验代码：
%欧拉法
function[x,y]=naeuler（dyfun,xspan,y0,h）
%dyfun是常微分方程，xspan是x的取值围，y0是初值，h是步长
x=xspan
（1）:
h:
xspan
（2）;
y
（1）=y0;
forn=1:
length（x）-1
y（n+1）=y（n）+h*feval（dyfun,x（n）,y（n））;
end
%改进的欧拉法
function[x,m,y]=naeuler2（dyfun,xspan,y0,h）
%dyfun是常微分方程，xspan是x的取值围，y0是初值，h是步长。
%返回值x为x取值，m为预报解，y为校正解
x=xspan
（1）:
h:
xspan
（2）;
y
（1）=y0;
m=zeros（length（x）-1,1）;
forn=1:
length（x）-1
k1=feval（dyfun,x（n）,y（n））;
y（n+1）=y（n）+h*k1;
m（n）=y（n+1）;
k2=feval（dyfun,x（n+1）,y（n+1））;
y（n+1）=y（n）+h*（k1+k2）/2;
end
%四阶S—K法
function[x,y]=rk（dyfun,xspan,y0,h）
%dyfun是常微分方程，xspan是x的取值围，y0是初值，h是步长。
x=xspan
（1）:
h:
xspan
（2）;
y
（1）=y0;
forn=1:
length（x）-1
k1=feval（dyfun,x（n）,y（n））;
k2=feval（dyfun,x（n）+h/2,y（n）+（h*k1）/2）;
k3=feval（dyfun,x（n）+h/2,y（n）+（h*k2）/2）;
k4=feval（dyfun,x（n）+h,y（n）+h*k3）;
y（n+1）=y（n）+（h/6）*（k1+2*k2+2*k3+k4）;
end
%主程序
x=[0:
0.1:
1];y=exp（-x）+x;
dyfun=inline（'-y+x+1'）;
[x1,y1]=naeuler（dyfun,[0,1],1,0.1）;
[x2,m,y2]=naeuler2（dyfun,[0,1],1,0.1）;
[x3,y3]=rk（dyfun,[0,1],1,0.1）;
plot（x,y,'r',x1,y1,'+',x2,y2,'*',x3,y3,'o'）;
xlabel（'x'）;ylabel（'y'）;
legend（'y为真解','y1为欧拉解','y2为改进欧拉解','y3为S—K解','Location','NorthWest'）;
1.2实验结果：
x
真解y
欧拉解y1
预报值m
校正值y2
S—K解y3
0.0
1.0000
1.0000
1.0000
1.0000
0.1
1.0048
1.0000
1.0000
1.0050
1.0048
0.2
1.0187
1.0100
1.0145
1.0190
1.0187
0.3
1.0408
1.0290
1.0371
1.0412
1.0408
0.4
1.0703
1.0561
1.0671
1.0708
1.0703
0.5
1.1065
1.0905
1.1037
1.1071
1.1065
0.6
1.1488
1.1314
1.1464
1.1494
1.1488
0.7
1.1966
1.1783
1.1945
1.1972
1.1966
0.8
1.2493
1.2305
1.2475
1.2500
1.2493
0.9
1.3066
1.2874
1.3050
1.3072
1.3066
1.0
1.3679
1.3487
1.3665
1.3685
1.3679
2、选取一种理论上收敛但是不稳定的算法对问题1进行计算，并与真解作比较。
（选改进的欧拉法）
2.1实验思路：
算法的稳定性是与步长h密切相关的。
而对于问题一而言，取定步长h=0.1不论是单步法或低阶多步法都是稳定的算法。
所以考虑改变h取值围，借此分析不同步长会对结果造成什么影响。
故依次采用h=2.0、2.2、2.4、2.6的改进欧拉法。
2.2实验代码：
%%主程序
x=[0:
3:
30];y=exp（-x）+x;
dyfun=inline（'-y+x+1'）;
[x1,m1,y1]=naeuler2（dyfun,[0,20],1,2）;
[x2,m2,y2]=naeuler2（dyfun,[0,22],1,2.2）;
[x3,m3,y3]=naeuler2（dyfun,[0,24],1,2.4）;
[x4,m4,y4]=naeuler2（dyfun,[0,26],1,2.6）;
subplot（2,2,1）
plot（x,y,'r',x1,y1,'+'）;xlabel（'h=2.0'）;
subplot（2,2,2）
plot（x,y,'r',x2,y2,'+'）;xlabel（'h=2.2'）;
subplot（2,2,3）
plot（x,y,'r',x3,y3,'+'）;xlabel（'h=2.4'）;
subplot（2,2,4）
plot（x,y,'r',x4,y4,'+'）;xlabel（'h=2.6'）;
2.3实验结果：
x
h=2.0
h=2.2
h=2.4
h=2.6
0.0
1.0000
1.0000
1.0000
1.0000
0.1
3.0000
3.4200
3.8800
4.3800
0.2
5.0000
5.8884
6.9904
8.3684
0.3
7.0000
8.4158
10.4418
13.4398
0.4
9.0000
11.0153
14.3979
20.4388
0.5
11.0000
13.7027
19.1008
30.8690
0.6
13.0000
16.4973
24.9092
47.4068
0.7
15.0000
19.4227
32.3536
74.8161
0.8
17.0000
22.5077
42.2194
121.5767
0.9
19.0000
25.7874
55.6687
202.7825
1.0
21.0000
29.3046
74.4217
345.3008
实验结果分析：
从实验1结果可以看出，在算法满足收敛性和稳定性的前提下，Eluer法虽然计算并不复杂，凡是精度不足，反观改进的Eluer法和S—K法虽然计算略微复杂但是结果很精确。
实验2改变了步长，导致算法理论上收敛但是不满足稳定性。
结果表示步长h越大，结果越失真。
对于同一个问题，步长h的选取变得尤为重要，这三种单步法的绝对稳定区间并不一样，所以并没有一种方法是万能的，我们应该根据不同的步长来选取合适的方法。
实验二：
Ritz-Galerkin方法与有限差分法
1、用中心差分格式求解边值问题
取步长h=0.1，并与真解作比较。
1.1实验代码：
%中心差分法
functionU=fdm（xspan,y0,y1,h）
%xspan为x取值围，y0,y1为边界条件，h为步长
N=1/h;
d=zeros（1,N-1）;
fori=1:
N
x（i）=xspan
（1）+i*h;
q（i）=1;
f（i）=x（i）;
end
fori=1:
N-1
d（i）=q（i）*h*h+2;
end
a=diag（d）;
b=zeros（N-1）;
c=zeros（N-1）;
fori=1:
N-2
b（i+1,i）=-1;
end
fori=1:
N-2
c（i,i+1）=-1;
end
A=a+b+c;
fori=2:
N-2
B（i,1）=f（i）*h*h;
end
B（1,1）=f
（1）*h*h+y0;
B（N-1,1）=f（N-1）*h*h+y1;
U=inv（A）*B;
%主程序
x=0:
0.1:
1;
y=x+（exp
（1）*exp（-x））/（exp
（2）-1）-（exp
（1）*exp（x））/（exp
（2）-1）;
y1=fdm（[0,1],0,0,0.1）;
y1=[0,y1',0];
plot（x,y,'r',x,y1,'+'）
xlabel（'x'）;ylabel（'y'）;
legend（'y真解','y1中心差分法','Location','NorthWest'）;
1.2实验结果：
x
y真解
y1中心差分法
0.0
0.0000
0.0000
0.1
0.0148
0.0148
0.2
0.0287
0.0287
0.3
0.0409
0.0408
0.4
0.0505
0.0504
0.5
0.0566
0.0565
0.6
0.0583
0.0582
0.7
0.0545
0.0545
0.8
0.0443
0.0443
0.9
0.0265
0.0265
1.0
0.0000
0.0000
2、用Ritz-Galerkin方法求解上述问题，并与真值作比较，列表画图。
2.1实验代码：
%Ritz_Galerkin法
functionvu=Ritz_Galerkin（x0,y0,x1,y1,h）
%x0,x1为x取值围，y0,y1为边界条件，h为步长
N=1/h;
symsx;
fori=1:
N
fai（i）=x*（1-x）*（x^（i-1））;
dfai（i）=diff（x*（1-x）*（x^（i-1）））;
end
fori=1:
N
forj=1:
N
fun=dfai（i）*dfai（j）+fai（i）*fai（j）;
A（i,j）=int（fun,x,0,1）;
end
fun=x*fai（i）+dfai（i）;
f（i）=int（fun,x,0,1）;
end
c=inv（A）*f';
product=c.*fai';
sum=0;
fori=1:
N
sum=sum+product（i）;
end
vu=[];
fory=0:
h:
1
v=subs（sum,x,y）;
vu=[vu,v];
end
y=0:
h:
1;
yy=0:
0.1:
1;
u=sin（yy）/sin
（1）-yy;
u=vpa（u,5）;
vu=vpa（vu,5）;
%主程序
x=0:
0.1:
1;
y=x+（exp
（1）*exp（-x））/（exp
（2）-1）-（exp
（1）*exp（x））/（exp
（2）-1）;
y1=Ritz_Galerkin（0,0,1,0,0.1）;
y1=double（y1）;
plot（x,y,'r',x,y1,'+'）
xlabel（'x'）;ylabel（'y'）;
legend（'y为真解','y1为R—G法','Location','NorthWest'）;
2.2实验结果：
x
y真解
y1R—G法
0.0
0.0000
0.0000
0.1
0.0148
0.0148
0.2
0.0287
0.0287
0.3
0.0409
0.0409
0.4
0.0505
0.0505
0.5
0.0566
0.0566
0.6
0.0583
0.0583
0.7
0.0545
0.0545
0.8
0.0443
0.0443
0.9
0.0265
0.0265
1.0
0.0000
0.0000
3、若边值条件为y（0）=0,y
（1）=1;则上述问题的数值解法怎样变化？
结果如何？
程序运算出来真解与数值解完全一样。
其值为y=x。
（具体运算不再赘述）。
实验结果分析：
对于实验1、2，我们可以看出不论是有限差分法还是Ritz-Galerkin法都非常稳定，结果非常精确（误差小于0.0001）。
对于实验3，编程中确定系数矩阵和常数项是最重要的。
确定过程中，要注意matlab中循环是从1开始的，而我们推导的公式中循环是从0开始的。
所以要区分开来谨慎对待，不然会产生极误差。
升级会员 