文档介绍：该【用Euler法改进的Euler法及四阶的龙格库塔法求解初值问题 】是由【莫比乌斯】上传分享，文档一共【10】页，该文档可以免费在线阅读，需要了解更多关于【用Euler法改进的Euler法及四阶的龙格库塔法求解初值问题 】的内容，可以使用淘豆网的站内搜索功能，选择自己适合的文档，以下文字是截取该文章内的部分文字，如需要获得完整电子版，请下载此文档到您的设备，方便您编辑和打印。微分方程数值解课程设计题目
1(30分)分别用Euler法、改进的Euler法及四阶的龙格库塔法求解初值问题:
,要求计算并绘出图形,最后比较三种算法的精度。
(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]
.当h=
>>Euler('fun',0,1,,10)
T=Columns1through10
Column11
②.当h=
>>Euler('fun',0,1,,5)
T=
(3).改进的Euler法:
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]
.当h=
>>Euler_modify('fun',0,1,,10)
T=Columns1through9
Columns10through11
②.当h=
>>Euler_modify('fun',0,1,,5)
T=
(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
.当h=
>>[x,y]=Rk_N4('fun',0,1,,10)
x=Columns1through9
Columns10through11
y=Columns1through9
Columns10through11
②.当h=
>>[x,y]=Rk_N4('fun',0,1,,5)
x=
y=
(5)作图:
1):h=
>>x=0::1;
>>y_euler=[];
>>y_euler_modify=[];
>>y_RK_N4=[];
plot(x,y_euler,'bo:',x,y_euler_modify,'r*-',x,y_RK_N4,'gv--');title('误差分析');xlabel('x轴');ylabel('y轴');text(,,'Euler法');text(,,'Euler改进法');text(,,'R-k法');text(,,'作者:李靖');text(,,'日期:');text(,,'h=');legend('Euler','Euler改进法','R_K法');gridon
2):h=
>>x=0::1;
>>y_euler=[];
>>y_euler_modify=[];
>>y_RK_N4=[];
plot(x,y_euler,'bo:',x,y_euler_modify,'r*-',x,y_RK_N4,'gv--');title('误差分析');xlabel('x轴');ylabel('y轴');text(,,'Euler法');text(,,'Euler改进法');text(,,'R-k法');text(,,'作者:xx');text(,,'日期:');text(,,'h=');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文件:
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文件互相的调用中,符号数组难以传值,,即结果为各点的导数的函数值。
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).:
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=
y=-
初始启动值为u(t1)=;u(t2)=
(2).Adams三步外插公式:
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,,,1/16,16)
x=Columns1through6
Columns7through12
Columns13through17
y=Columns1through6
--
Columns7through12
------
Columns13through17
-----
(3).Adams四步外插公式:
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,,,,1/16,16)
x=Columns1through6
Columns7through12
Columns13through17
y=Columns1through6
--
Columns7through12
------
Columns13through17
-----
(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
--
Columns7through12
------
Columns13through17
-----
(5).作图分析误差:
>>x=0:1/16:1;
>>yadams3=[-------------];
>>yadams4=[-------------];
>>y_exact=[-------------];
plot(x,yadams3,'bo:',x,yadams4,'r*-',x,y_exact,'k-');title('误差分析');text(,0,'Adams三步公式');text(,-,'Adams四步公式');text(,-,'精确解');text(,,'作者:李靖');text(,,'日期:');legend('Adams三步公式','Adams四步公式','精确解');gridon
从图中我们可以清楚的看到三步公式和四步公式在[0,]之间与精确解的值相差无几,并且四步公式比三步公式精确度更高。但是在[,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
Columns15through28
Columns29through42
Columns43through56
Columns57through70
Columns71through81