计算方法上机作业——龙格库塔法matlab程序.pdf

《计算方法上机作业——龙格库塔法matlab程序.pdf》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、附录5龙格—库塔法求解常微分方程初值问题的matlab程序附录5龙格—库塔法求解常微分方程初值问题的matlab程序clear;clc;%======需要输入的=======%n=input('请输入一阶微分方程组的个数或者高阶方程组的阶数n：');a=input('请输入求解区间的下限a:');b=input('请输入求解区间的上限b:');h=input('请输入求解步长h:');fprintf('请依次输入%.f个方程的右端函数fi(x,y(x)):'

2、,n)s=cell(n,1);forp=1:nfprintf('f%.f(x,y(x))=',p)r=input('','s');s{p}=r;endfprintf('请依次输入%.f个初值yi(%.f):',n,a)y0=zeros(n,1);forq=1:nfprintf('y%.f(%.f)=',q,a);y0(q)=input('');end%======需要输入的=======%N=(b-a)/h+1;x=a:h:b;K=zeros(n,4);y=zeros(n,N);y(:,1

3、)=y0;fori=1:N-1fora=1:nf=inline(char(s(a)),'x','y');K(a,1)=h*f(x(i),y(:,i));endforb=1:nf=inline(char(s(b)),'x','y');K(b,2)=h*f(x(i)+h/2,y(:,i)+K(:,1)/2);endforc=1:nf=inline(char(s(c)),'x','y');K(c,3)=h*f(x(i)+h/2,y(:,i)+K(:,2)/2);endford=1:nf=inline(c

4、har(s(d)),'x','y');K(d,4)=h*f(x(i)+h,y(:,i)+K(:,3));endforj=1:ny(j,i+1)=y(j,i)+1/6*(K(j,1)+2*K(j,2)+2*K(j,3)+K(j,4));endplot(x(i),y(1,i),'.k','markersize',15)holdon;endgrid;fprintf('%10c','x')28计算方法上机报告formm=1:nfprintf('%9c%.f','y',mm)endfornn=1:Nfp

5、rintf('%10.5f',x(nn));forll=1:nfprintf('%10.5f',y(ll,nn));endifnn==Nfprintf('');endend29