用Euler法改进的Euler法及四阶的龙格库塔法求解初值问题.docx
- 文档编号:6324750
- 上传时间:2023-01-05
- 格式:DOCX
- 页数:12
- 大小:45.96KB
用Euler法改进的Euler法及四阶的龙格库塔法求解初值问题.docx
《用Euler法改进的Euler法及四阶的龙格库塔法求解初值问题.docx》由会员分享,可在线阅读,更多相关《用Euler法改进的Euler法及四阶的龙格库塔法求解初值问题.docx(12页珍藏版)》请在冰豆网上搜索。
用Euler法改进的Euler法及四阶的龙格库塔法求解初值问题
微分方程数值解课程设计题目
1(30分)分别用Euler法、改进的Euler法及四阶的龙格库塔法求解初值问题:
h分别取0.1和0.2,要求计算并绘出图形,最后比较三种算法的精度。
(1).先构建初值问题的函数
M文件:
functionz=fun(x,y)
z=y-2*x/y;
(2).Euler法
M文件:
functionE=Euler(fun,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x
(1)=x0;y
(1)=y0;
forn=1:
N
x(n+1)=x(n)+h;
y(n+1)=y(n)+h*feval(fun,x(n),y(n));
end
T=[x;y]
1.当h=0.1时
>>Euler('fun',0,1,0.1,10)
T=Columns1through10
00.10000.20000.30000.40000.50000.60000.70000.80000.9000
1.00001.10001.19181.27741.35821.43511.50901.58031.64981.7178
Column11
1.0000
1.7848
②.当h=0.2时
>>Euler('fun',0,1,0.2,5)
T=00.20000.40000.60000.80001.0000
1.00001.20001.37331.53151.68111.8269
(3).改进的Euler法:
M文件:
Euler_modify.m
functionE=Euler_modify(fun,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x
(1)=x0;y
(1)=y0;
forn=1:
N
x(n+1)=x(n)+h;
z0=y(n)+h*feval(fun,x(n),y(n));
y(n+1)=y(n)+h/2*(feval(fun,x(n),y(n))+feval(fun,x(n+1),z0));
end
T=[x;y]
1.当h=0.1时
>>Euler_modify('fun',0,1,0.1,10)
T=Columns1through9
00.10000.20000.30000.40000.50000.60000.70000.8000
1.00001.09591.18411.26621.34341.41641.48601.55251.6165
Columns10through11
0.90001.0000
1.67821.7379
②.当h=0.2时
>>Euler_modify('fun',0,1,0.2,5)
T=00.20000.40000.60000.80001.0000
1.00001.18671.34831.49371.62791.7542
(4).四阶R_K方法:
M文件:
function[x,y]=Rk_N4(f,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x
(1)=x0;y
(1)=y0;
forn=1:
N
x(n+1)=x(n)+h;
k1=h*feval(f,x(n),y(n));
k2=h*feval(f,x(n)+1/2*h,y(n)+1/2*k1);
k3=h*feval(f,x(n)+1/2*h,y(n)+1/2*k2);
k4=h*feval(f,x(n)+h,y(n)+k3);
y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);
end
1.当h=0.1时
>>[x,y]=Rk_N4('fun',0,1,0.1,10)
x=Columns1through9
00.10000.20000.30000.40000.50000.60000.70000.8000
Columns10through11
0.90001.0000
y=Columns1through9
1.00001.09541.18321.26491.34161.41421.48321.54921.6125
Columns10through11
1.67331.7321
②.当h=0.2时
>>[x,y]=Rk_N4('fun',0,1,0.2,5)
x=00.20000.40000.60000.80001.0000
y=1.00001.18321.34171.48331.61251.7321
(5)作图:
1):
h=0.1
>>x=0:
0.1:
1;
>>y_euler=[1.00001.10001.19181.27741.35821.43511.50901.58031.64981.71781.7848];
>>y_euler_modify=[1.00001.09591.18411.26621.34341.41641.48601.55251.61651.67821.7379];
>>y_RK_N4=[1.00001.09541.18321.26491.34161.41421.48321.54921.61251.67331.7321];
plot(x,y_euler,'bo:
',x,y_euler_modify,'r*-',x,y_RK_N4,'gv--');title('误差分析');xlabel('x轴');ylabel('y轴');text(0.3,1.3,'Euler法');text(0.4,1.35,'Euler改进法');text(0.5,1.4,'R-k法');text(0.8,1.02,'作者:
李靖');text(0.8,1.07,'日期:
2010.6.19');text(0.4,1.75,'h=0.1');legend('Euler','Euler改进法','R_K法');gridon
2):
h=0.2
>>x=0:
0.2:
1;
>>y_euler=[1.00001.20001.37331.53151.68111.8269];
>>y_euler_modify=[1.00001.18671.34831.49371.62791.7542];
>>y_RK_N4=[1.00001.18321.34171.48331.61251.7321];
plot(x,y_euler,'bo:
',x,y_euler_modify,'r*-',x,y_RK_N4,'gv--');title('误差分析');xlabel('x轴');ylabel('y轴');text(0.3,1.3,'Euler法');text(0.4,1.35,'Euler改进法');text(0.5,1.4,'R-k法');text(0.8,1.02,'作者:
xx');text(0.8,1.07,'日期:
2010.6.19');text(0.4,1.75,'h=0.2');legend('Euler','Euler改进法','R_K法');gridon
2.(20分)编写一个程序用tylor级数法求解问题:
取tylor级数法的截断误差为O(
),即要用u(t),u’(t),…,u(20)(t)的值。
(提示:
可用一个简单的递推公式来获得u(n)(t),n=1,2,……)。
(1).先利用递推公式计算u(t)的各阶导数
M文件:
differ.m
functionD=differ(N)
y=zeros(1,N+1);
symsxy;
y
(1)=x*y;
forn=1:
N
y(n+1)=y(n)*x+diff(y(n),x);
end
y
输入(由于随着导数的阶数增高,项越来越多。
所以只举例计算到5阶):
>>differ(5)
y=
[x*y,x^2*y+y,(x^2*y+y)*x+2*x*y,((x^2*y+y)*x+2*x*y)*x+3*x^2*y+3*y,(((x^2*y+y)*x+2*x*y)*x+3*x^2*y+3*y)*x+(3*x^2*y+3*y)*x+(x^2*y+y)*x+8*x*y,((((x^2*y+y)*x+2*x*y)*x+3*x^2*y+3*y)*x+(3*x^2*y+3*y)*x+(x^2*y+y)*x+8*x*y)*x+((3*x^2*y+3*y)*x+(x^2*y+y)*x+8*x*y)*x+((x^2*y+y)*x+2*x*y)*x+15*x^2*y+15*y]
(2).利用tylor级数法计算初值问题:
由于M文件互相的调用中,符号数组难以传值,所以我就将上面的differ.m改为赋值的计算,即结果为各点的导数的函数值。
M文件1:
functionD=tylor_differ(x0,y0,N)
symsxy;
t1=x*y;
t2=x*y;
forn=1:
N-1
t2=t1*x+diff(t1,x);
t1=t2;
end
D=subs(t2,{x,y},{x0,y0});
M文件2:
迭代计算微分方程初值问题
functionT=Tylor(x0,y0,h,N)
t=zeros(1,N+1);u=zeros(1,N+1);
t
(1)=x0;u
(1)=y0;
symsxy;
forn=1:
N
t(n+1)=t(n)+h;w=1;
fors=1:
20%这里可设置取tylor级数的多少项,即要函数的倒数的阶数
w=w*s;
tempt=h^s/w*tylor_differ(t(n),u(n),s);
u(n+1)=u(n+1)+tempt;
end
u(n+1)=u(n+1)+u(n);
end
T=[t;u]
输入:
3(30分)分别用Adams三步和四步外插公式,h=1/16求解
将计算结果与真解u(t)=e-8t+t2/2-t进行比较,并对所产生的现象进行理论分析。
(1).构建初值函数的M文件fa.m:
functionz=fa(x,y)
z=-8*y+4*x^2-7*x-1;
先利用R_K四阶方法计算Adams公式中的插值点fn,fn-1,fn-2即初始值u(t0),u(t1),u(t2)来作为启动值
>>[x,y]=Rk_N4('fa',0,1,1/16,4)
x=00.06250.12500.18750.2500
y=1.00000.54620.25100.0535-0.0832
初始启动值为u(t1)=0.5462;u(t2)=0.2510
(2).Adams三步外插公式:
M文件:
Adams3.m
function[x,y]=Adams3(fun,x0,y0,y1,y2,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x
(1)=x0;y
(1)=y0;
x
(2)=x0+h;y
(2)=y1;
x(3)=x0+2*h;y(3)=y2;
forn=3:
N
x(n+1)=x(n)+h;
f1=feval(fun,x(n),y(n));
f2=feval(fun,x(n-1),y(n-1));
f3=feval(fun,x(n-2),y(n-2));
y(n+1)=y(n)+h/12*(23*f1-16*f2+5*f3);
end
>>[x,y]=Adams3('fa',0,1,0.5462,0.2510,1/16,16)
x=Columns1through6
00.06250.12500.18750.25000.3125
Columns7through12
0.37500.43750.50000.56250.62500.6875
Columns13through17
0.75000.81250.87500.93751.0000
y=Columns1through6
1.00000.54620.25100.0416-0.0909-0.1940
Columns7through12
-0.2606-0.3202-0.3592-0.3984-0.4234-0.4503
Columns13through17
-0.4658-0.4830-0.4904-0.4990-0.4987
(3).Adams四步外插公式:
M文件:
Adams4.m
function[x,y]=Adams4(fun,x0,y0,y1,y2,y3,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x
(1)=x0;y
(1)=y0;
x
(2)=x0+h;y
(2)=y1;
x(3)=x0+2*h;y(3)=y2;
x(4)=x0+3*h;y(4)=y3;
forn=4:
N
x(n+1)=x(n)+h;
f1=feval(fun,x(n),y(n));
f2=feval(fun,x(n-1),y(n-1));
f3=feval(fun,x(n-2),y(n-2));
f4=feval(fun,x(n-3),y(n-3));
y(n+1)=y(n)+h/24*(55*f1-59*f2+37*f3-9*f4);
end
>>[x,y]=Adams4('fa',0,1,0.5462,0.2510,0.0535,1/16,16)
x=Columns1through6
00.06250.12500.18750.25000.3125
Columns7through12
0.37500.43750.50000.56250.62500.6875
Columns13through17
0.75000.81250.87500.93751.0000
y=Columns1through6
1.00000.54620.25100.0535-0.0790-0.1795
Columns7through12
-0.2483-0.3124-0.3487-0.3996-0.4102-0.4630
Columns13through17
-0.4417-0.5151-0.4414-0.5687-0.3970
(4).精确解:
M文件:
functionE=fa_analytic(x0,y0,h,N)
symsx;
f=1/2*x^2-x+exp(-8*x);
temp=zeros(1,N+1);
forn=1:
N+1
temp(n)=subs(f,x,x0);
x0=x0+h;
end
T=temp
>>fa_analytic(0,1,1/16,16)
T=Columns1through6
1.00000.54600.25070.0532-0.0834-0.1816
Columns7through12
-0.2549-0.3116-0.3567-0.3932-0.4229-0.4471
Columns13through17
-0.4663-0.4809-0.4913-0.4975-0.4997
(5).作图分析误差:
>>x=0:
1/16:
1;
>>yadams3=[1.00000.54620.25100.0416-0.0909-0.1940-0.2606-0.3202-0.3592-0.3984-0.4234-0.4503-0.4658-0.4830-0.4904-0.4990-0.4987];
>>yadams4=[1.00000.54620.25100.0535-0.0790-0.1795-0.2483-0.3124-0.3487-0.3996-0.4102-0.4630-0.4417-0.5151-0.4414-0.5687-0.3970];
>>y_exact=[1.00000.54600.25070.0532-0.0834-0.1816-0.2549-0.3116-0.3567-0.3932-0.4229-0.4471-0.4663-0.4809-0.4913-0.4975-0.4997];
plot(x,yadams3,'bo:
',x,yadams4,'r*-',x,y_exact,'k-');title('误差分析');text(0.2,0,'Adams三步公式');text(0.3,-0.17,'Adams四步公式');text(0.5,-0.37,'精确解');text(0.7,0.8,'作者:
李靖');text(0.7,0.7,'日期:
2010.6.19');legend('Adams三步公式','Adams四步公式','精确解');gridon
从图中我们可以清楚的看到三步公式和四步公式在[0,0.5]之间与精确解的值相差无几,并且四步公式比三步公式精确度更高。
但是在[0.5,1]之间时四步公式出现了波动。
4(20分)选择一个精度不低于四阶的方法,对t>=0时的标准正态分布函数:
产生一张在[0,5]之间的80个等距节点(即h=1/16)处的函数值表。
提示寻找一个以Φ(x)为解的初值问题。
(1).先建立微分方程初值问题:
M文件:
functionz=normsdist(x,y)
z=1/sqrt(2*pi)*exp(-x/2);
(2).利用四阶R_K方法计算微分方程初值问题:
>>[x,y]=Rk_N4('normsdist',0,1,1/16,80)
y=Columns1through14
1.00001.02451.04831.07141.09381.11541.13641.15681.17651.19561.21411.23211.24951.2664
Columns15through28
1.28271.29861.31391.32881.34331.35731.37081.38391.39671.40901.42101.43261.44381.4547
Columns29through42
1.46531.47551.48541.49501.50441.51341.52211.53061.53881.54681.55451.56201.56931.5763
Columns43through56
1.58311.58971.59611.60241.60841.61421.61991.62531.63061.63581.64081.64561.65031.6548
Columns57through70
1.65921.66351.66761.67161.67551.67931.68291.68651.68991.69321.69641.69961.70261.7055
Columns71through81
1.70841.71111.71381.71641.71891.72131.72371.72601.72821.73031.7324
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Euler 改进 龙格库塔法 求解 初值问题