四阶龙格-库塔法求解常微分方程的初值问题-matlab通用程序

参考教材《数值分析》李乃成.梅立泉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('请输入在初值点的函数值f(a)=');x=a;y11(1)=y(1,1);fork1=2:(mm+1)y1=y(1,1);K(1,1)=h*(eval(ym));%计算K1x=x+h/2;y(1,1)=y1+K(1,1)/2;y1=y(1,1);K(1,2)=h*(eval(ym));%计算K2x=x;y(1,1)=y1+K(1,2)/2-K(1,1)/2;y1=y(1,1);K(1,3)=h*(eval(ym));%计算K3x=x+h/2;y(1,1)=y1+K(1,3)-K(1,2)/2;y1=y(1,1);K(1,4)=h*(eval(ym));%计算K4y11(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;%一共要求解mm个数据点fork2=1:m%读取初值条件fprintf('请输入%d阶导数的初值f(%d)(a)=\n',(k2-1),(k2-1));y(k2,1)=input('=');endfork2=1:my22(1,k2)=y(k2,1);%先把初值保存在矩阵y22(m,n)中,m表示第几个所求点，n表示第n阶初值endx=a;fork4=2:(mm+1)%求解mm个数据点的循环fork=1:(m-1)%计算K1，包括每一阶的K1K(k,1)=h*y(k+1,1);%y(k+1,1)中k+1表示第k+1阶，1表示第一个点；K(k,1)中k表示阶数，1表示K1endK(m,1)=h*(eval(ym));x=x+h/2;%求解K1之前，先重新对x和y赋值fork3=1:my(k3,1)=y(k3,1)+K(k3,1)/2;endfork=1:(m-1)%计算K2K(k,2)=h*y(k+1,1);endK(m,2)=h*(eval(ym));x=x;fork3=1:my(k3,1)=y(k3,1)-K(k3,1)/2+K(k3,2)/2;endfork=1:(m-1)%计算K3K(k,3)=h*y(k+1,1);endK(m,3)=h*(eval(ym));x=x+h/2;fork3=1:my(k3,1)=y(k3,1)+K(k3,3)-K(k3,2)/2;%这里容易出错endfork=1:(m-1)%计算K4K(k,4)=h*y(k+1,1);endK(m,4)=h*(eval(ym));fork5=1:my22(k4,k5)=y22(k4-1,k5)+(K(k5,1)+2*K(k5,2)+2*K(k5,3)+K(k5,4))/6;%这里，除了要求出下一个点的数值，还要求出相应的导数值endfork6=1:m%除了对y(1,1)重新赋值外，还要对y(2,1)等重新赋值y(k6,1)=y22(k4,k6);endx=a+(k4-1)*h;endy22(:,1)end