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

精品文档四阶龙格库塔法解微分方程一、四阶龙格库塔法解一阶微分方程一阶微分方程： ，初始值y(0)=0，求解区间0 10。MATLAB程序：% 四阶龙哥库塔法解一阶微分方程% y=cost% y(0)=0, 0t10，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); for i=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;end subplot(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一阶微分方程： ，初始值x(1)=2，求解区间1 3。MATLAB程序：% 四阶龙哥库塔法解微分方程% x(t)=(t*x-x2)/t2% x(1)=2, 1t3, 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; for i=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;end subplot(211)plot(t,x,-);xlabel(t);ylabel(x);legend(Approximation); % ode45求精确解t0=t(1);x0=x(1);xspan=t0 tf;x_ode45,y_ode45=ode45(F_tx,xspan,x0);subplot(212)plot(x_ode45,y_ode45,-); xlabel(t);ylabel(x);legend(Exact);图2二、四阶龙格库塔法解二阶微分方程二阶微分方程： ，初始值y(0)=0，y(0)=-1，求解区间0 10。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); for i=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);end plot(t,y)xlabel(t);ylabel(y);title(y=sin(t);legend(y=sint);图3二阶微分方程： ，初始值y(1)=0，y(1)=0，求解区间0 25。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; for i=1:length(x)-1 K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i); % K1 and L1 K2=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); % K4 and L4 RK_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);end plot(x,RK_y,b-);xlabel(Variable x);ylabel(Variable y);A=x,RK_y;y,T=max(RK_y);legend(RK4方法);fprintf(角度峰值等于%d,y) %角度的峰值也就是fprintf(n)fprintf(周期等于%d,T*h) %周期function w=rightf_sys1(x,y