四阶龙格_库塔法求解常微分方程的初值问题_matlab通用程序.docxVIP

参考教材《数值分析》李乃成.梅立泉 clearclcformat longm=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);if m==1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? %一阶初值问题单独求解? ? mm=(b-a)/h;? ? y(1,1)=input(请输入在初值点的函数值f(a)=);? ? x=a;? ? y11(1)=y(1,1);? ? for k1=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)); ? ? ? ? ? ? ? ? ? ? ? ? ?%计算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; ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? %一共要求解mm个数据点? ? for k2=1:m ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?%读取初值条件? ? ? ? fprintf(请输入%d阶导数的初值f(%d)(a)=\n,(k2-1),(k2-1));? ? ? ? y(k2,1)=input(=);? ? end? ? for k2=1:m ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ?y22(1,k2)=y(k2,1); ? ? ? ? ? ? ? ? ? ? ? ? ?%先把初值保存在矩阵y22(m,n)中,m表示第几个所求点，n表示第n阶初值? ? end? ? x=a;? ? for k4=2:(mm+1) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? %求解mm个数据点的循环? ? ? ? for k=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赋值? ? ? ? for k3=1:m ? ? ? ? ? ? ?? ? ? ? ? ? y(k3,1)=y(k3,1)+K(k3,1)/2;? ? ? ? end? ? ? ? for k=1:(m-1) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?%计算K2? ? ? ? ? ? K(k,2)=h*y(k+1,1);? ? ? ? end? ? ? ? K(m,2)=h*(eval(ym));? ? ? ? x=x;? ? ? ? for k3=1:m