欧拉法RK法的MATLAB程序微分方程数值解
dy/dx=y/2x+x2/2y
y(1)=1
a) 欧拉格式
代码:
function f=ww(x,y)
f=y/(2*x)+x^2/(2*y);
h=0.1; y(1)=1; x(1)=1;
n=1;
while x(n)<=2
n=n+1;
x(n)=x(n-1)+h;
y(n)=y(n-1)+h*ww(x(n-1),y(n-1));
x(n)
y(n)
end
显示结果: xn yn
1...
微分方程数值解
dy/dx=y/2x+x2/2y
y(1)=1
a) 欧拉格式
代码:
function f=ww(x,y)
f=y/(2*x)+x^2/(2*y);
h=0.1; y(1)=1; x(1)=1;
n=1;
while x(n)<=2
n=n+1;
x(n)=x(n-1)+h;
y(n)=y(n-1)+h*ww(x(n-1),y(n-1));
x(n)
y(n)
end
显示结果: xn yn
1.0000 1.0000
1.1000 1.1000
1.2000 1.2050
1.3000 1.3150
1.4000 1.4298
1.5000 1.5494
1.6000 1.6737
1.7000 1.8024
1.8000 1.9356
1.9000 2.0731
2.0000 2.2147
若h=0.05 y=2.2254 ;h=0.01 y=2.2339
b)中点格式
代码:
h=0.1; y(1)=1; x(1)=1;
n=1;
while x(n)<=2
n=n+1;
x(n)=x(n-1)+h;
y(n)=y(n-1)+h*ww(x(n-1)+h/2,y(n-1)+h*ww(x(n-1),y(n-1))/2);
x(n)
y(n)
end
显示结果xn yn
1.0000 1.0000
1.1000 1.1025
1.2000 1.2100
1.3000 1.3223
1.4000 1.4395
1.5000 1.5613
1.6000 1.6877
1.7000 1.8184
1.8000 1.9535
1.9000 2.0928
2.0000 2.2362
c)预报校正格式
代码:
h=0.1;y(1)=1;x(1)=1;
n=1;
while x(n)<=2
n=n+1;
x(n)=x(n-1)+h;
yy(n)=y(n-1)+h*ww(x(n-1),y(n-1));
y(n)=y(n-1)+h*[ww(x(n-1),y(n-1))+ww(x(n),yy(n))]/2;
x(n)
y(n)
end
显示结果:xn yn
1.0000 1.0000
1.1000 1.1025
1.2000 1.2100
1.3000 1.3223
1.4000 1.4395
1.5000 1.5613
1.6000 1.6876
1.7000 1.8184
1.8000 1.9535
1.9000 2.0928
2.0000 2.2361
d)R-K格式
代码:
h=0.1;y(1)=1;x(1)=1;
n=1;
while x(n)<=2
n=n+1;
x(n)=x(n-1)+h;
k(1)=ww(x(n-1),y(n-1));
k(2)=ww(x(n-1)+h/2,y(n-1)+h*k(1)/2);
k(3)=ww(x(n-1)+h/2,y(n-1)+h*k(2)/2);
k(4)=ww(x(n-1)+h/2,y(n-1)+h*k(3));
y(n)=y(n-1)+h*(k(1)+2*k(2)+2*k(3)+k(4))/6;
x(n)
y(n)
end
显示结果;xn yn
1.0000 1.0000
1.1000 1.1021
1.2000 1.2091
1.3000 1.3210
1.4000 1.4377
1.5000 1.5590
1.6000 1.6849
1.7000 1.8152
1.8000 1.9498
1.9000 2.0886
2.0000 2.2316
2008-05-15
本文档为【欧拉法RK法的MATLAB程序】,请使用软件OFFICE或WPS软件打开。作品中的文字与图均可以修改和编辑,
图片更改请在作品中右键图片并更换,文字修改请直接点击文字进行修改,也可以新增和删除文档中的内容。
该文档来自用户分享,如有侵权行为请发邮件ishare@vip.sina.com联系网站客服,我们会及时删除。
[版权声明] 本站所有资料为用户分享产生,若发现您的权利被侵害,请联系客服邮件isharekefu@iask.cn,我们尽快处理。
本作品所展示的图片、画像、字体、音乐的版权可能需版权方额外授权,请谨慎使用。
网站提供的党政主题相关内容(国旗、国徽、党徽..)目的在于配合国家政策宣传,仅限个人学习分享使用,禁止用于任何广告和商用目的。