资源描述：
《四阶龙格-库塔法求解常微分方程的初值问题-matlab通用程序.doc》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。
1、参考教材《数值分析》李乃成.梅立泉clearclcformatlongm=input('请输入常微分方程的阶数m=');a=input('请输入x下限a=');b=input('请输入x上限b=');h=input('请输入步长h=');ym=input('令y(1,1)=y,y（2,1）=y’,y(3,1)=y’’...请输入ym=','s'); %输入的时候必须按照这个形式输入y1=y(1,1);ifm==1 %一阶初值问题单独求解 mm=(b-a)/h; y(1,1)=input('请输入在初值点的函
2、数值f(a)='); x=a; y11(1)=y(1,1); fork1=2:(mm+1) y1=y(1,1); K(1,1)=h*(eval(ym)); %计算K1 x=x+h/2; y(1,1)=y1+K(1,1)/2; y1=y(1,1); K(1,2)=h*(eval(ym)); %计算K2 x=x; y(1,1)=y1+K(1,2)/2-K(1,1)/2; y1=y(1,1); K(1,3)=h*(eval(ym));
3、 %计算K3 x=x+h/2; y(1,1)=y1+K(1,3)-K(1,2)/2; y1=y(1,1); K(1,4)=h*(eval(ym)); %计算K4 y11(k1)=y11(k1-1)+(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4))/6; y(1,1)=y11(k1); x=a+(k1-1)*h; endy11else %高阶初值问题 mm=(b-a)/h; %一共
4、要求解mm个数据点 fork2=1:m %读取初值条件 fprintf('请输入%d阶导数的初值f(%d)(a)=',(k2-1),(k2-1)); y(k2,1)=input('='); end fork2=1:m y22(1,k2)=y(k2,1); %先把初值保存在矩阵y22(m,n)中,m表示第几个所求点，n表示第n阶初值 end x=a; fork4=2:(mm+1) %求解mm个
5、数据点的循环 fork=1:(m-1) %计算K1，包括每一阶的K1 K(k,1)=h*y(k+1,1); %y(k+1,1)中k+1表示第k+1阶，1表示第一个点；K(k,1)中k表示阶数，1表示K1 end K(m,1)=h*(eval(ym)); x=x+h/2; %求解K1之前，先重新对x和y赋值 fork3=1:m y(k3,1)=y(k3,1)+K(k3,1)/2; end fork
6、=1:(m-1) %计算K2 K(k,2)=h*y(k+1,1); end K(m,2)=h*(eval(ym)); x=x; fork3=1:m y(k3,1)=y(k3,1)-K(k3,1)/2+K(k3,2)/2; end fork=1:(m-1) %计算K3 K(k,3)=h*y(k+1,1); end K(m,3)=h*(eval(ym)); x=x+h/2; fork3=1:m y(
7、k3,1)=y(k3,1)+K(k3,3)-K(k3,2)/2; %这里容易出错 end fork=1:(m-1) %计算K4 K(k,4)=h*y(k+1,1); end K(m,4)=h*(eval(ym)); fork5=1:m y22(k4,k5)=y22(k4-1,k5)+(K(k5,1)+2*K(k5,2)+2*K(k5,3)+K(k5,4))/6; %这里，除了要求出下一个点的数值，还要求出相应的导数值 end fork6=1:m
8、 %除了对y(1,1)重新赋值外，还要对y
