四阶龙格库塔法解微分方程

四阶龙格库塔法解微分方程一、四阶龙格库塔法解一阶微分方程①一阶微分方程：cosyt，初始值y(0)=0，求解区间[010]。MATLAB程序：%%%%%%%%%%%四阶龙哥库塔法解一阶微分方程%%%%%%%%%%%y'=cost%%%%%%%%%%%y(0)=0,0≤t≤10，h=0.01%%%%%%%%%%%y=sinth=0.01;hf=10;t=0:h:hf;y=zeros(1,length(t));y(1)=0;F=@(t,y)(cos(t));fori=1:(length(t)-1)k1=F(t(i),y(i));k2=F(t(i)+h/2,y(i)+k1*h/2);k3=F(t(i)+h/2,y(i)+k2*h/2);k4=F(t(i)+h,y(i)+k3*h);y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4)*h;endsubplot(211)plot(t,y,'-')xlabel('t');ylabel('y');title('Approximation');span=[0,10];p=y(1);[t1,y1]=ode45(F,span,p);subplot(212)plot(t1,y1)xlabel('t');ylabel('y');title('Exact');图1②一阶微分方程：22*/xtxxt，初始值x(1)=2，求解区间[13]。MATLAB程序：%%%%%%%%%%%四阶龙哥库塔法解微分方程%%%%%%%%%%%x'(t)=(t*x-x^2)/t^2%%%%%%%%%%%x(1)=2,1≤t≤3,h=1/128%%%%%%%%%%%精确解：x(t)=t/(0.5+lnt)h=1/128;%%%%%步长tf=3;t=1:h:tf;x=zeros(1,length(t));x(1)=2;%%%%%初始值F_tx=@(t,x)(t.*x-x.^2)./t.^2;fori=1:(length(t)-1)k_1=F_tx(t(i),x(i));k_2=F_tx(t(i)+0.5*h,x(i)+0.5*h*k_1);k_3=F_tx((t(i)+0.5*h),(x(i)+0.5*h*k_2));k_4=F_tx((t(i)+h),(x(i)+k_3*h));x(i+1)=x(i)+(1/6)*(k_1+2*k_2+2*k_3+k_4)*h;endsubplot(211)plot(t,x,'-');xlabel('t');ylabel('x');legend('Approximation');%%%%%%%%%%%%%%%%%%%%%%%%%%%%ode45求精确解t0=t(1);x0=x(1);xspan=[t0tf];[x_ode45,y_ode45]=ode45(F_tx,xspan,x0);subplot(212)plot(x_ode45,y_ode45,'--');xlabel('t');ylabel('x');legend('Exact');图2二、四阶龙格库塔法解二阶微分方程①二阶微分方程：cosyt，初始值y(0)=0，y'(0)=-1，求解区间[010]。MATLAB程序：%%%%%%%%%龙格库塔欧拉方法解二阶微分方程%%%%%%%%%y''=cost%%%%%%%%%y'(0)=-1,y(0)=0%%%%%%%%%转换为一阶微分方程h=0.01;ht=10;t=0:h:ht;%%%%%%%%%y'=z1,z1'=y''=costy=zeros(1,length(t));z=zeros(1,length(t));y(1)=0;z(1)=-1;f=@(t,y,z)(z);g=@(t,y,z)(sin(t));fori=1:(length(t)-1)K1=f(t(i),y(i),z(i));L1=g(t(i),y(i),z(i));K2=f(t(i)+h/2,y(i)+h/2*K1,z(i)+h/2*L1);L2=g(t(i)+h/2,y(i)+h/2*K1,z(i)+h/2*L1);K3=f(t(i)+h/2,y(i)+h/2*K2,z(i)+h/2*L2);L3=g(t(i)+h/2,y(i)+h/2*K2,z(i)+h/2*L2);K4=f(t(i)+h,y(i)+h*K3,z(i)+h*L3);L4=g(t(i)+h,y(i)+h*K3,z(i)+h*L3);y(i+1)=y(i)+h/6*(K1+2*K2+2*K3+K4);z(i+1)=z(i)+h/6*(L1+2*L2+2*L3+L4);endplot(t,y)xlabel('t');ylabel('y');title('y''=sin(t)');legend('y''=sint');图3②二阶微分方程：7.35499*cosyx，初始值y(1)=0，y'(1)=0，求解区间[025]。MATLAB程序：%%%%%%%%%龙格库塔欧拉方法解二阶微分方程%%%%%%%%%y''=7.35499*cosx%set(0,'RecursionLimit',500)h=0.01;a=0;b=25;x=a:h:b;RK_y(1)=0;%初值RK_z(1)=0;fori=1:length(x)-1K1=RK_z(i);L1=rightf_sys1(x(i),RK_y(i),RK_z(i));%K1andL1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);K4=RK_z(i)+h*L3;L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);%K4andL4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,RK_y,'b-');xlabel('Variablex');ylabel('Variabley');A=[x,RK_y];[y,T]=max(RK_y);legend('RK4方法');fprintf('角度峰值等于%d',y)%角度的峰值也就是πfprintf('\n')fprintf('周期等于%d',T*h)%周期functionw=rightf_sys1(x,y,z)w=7.35499*cos(y);end图4注：四阶龙格库塔法求解二阶微分方程时，只需把二阶微分方程转化为一阶微分方程，再进行相关求解即可。