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

《四阶龙格库塔法解微分方程.docx》由会员分享，可在线阅读，更多相关《四阶龙格库塔法解微分方程.docx（8页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、四阶龙格库塔法解微分方程四阶龙格库塔法解微分方程一、四阶龙格库塔法解一阶微分方程一阶微分方程：cosyt&，初始值y(0)=0，求解区间010。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);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)+k

2、3*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-x2)/t2%x(1)=2,1t3,h=1/1

3、28%准确解：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

4、(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,-);%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

5、=(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(

6、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:len

7、gth(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);