四階龍格庫塔法解微分方程(共6頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上四階龍格庫塔法解微分方程一、四階龍格庫塔法解一階微分方程一階微分方程： ，初始值y(0)=0，求解區(qū)間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(

2、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一階微分方程： ，

3、初始值x(1)=2，求解區(qū)間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

4、)+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

5、_ode45,'-'); xlabel('t');ylabel('x');legend('Exact');圖2二、四階龍格庫塔法解二階微分方程二階微分方程： ，初始值y(0)=0，y'(0)=-1，求解區(qū)間0 10。MATLAB程序：% 龍格庫塔歐拉方法解二階微分方程% y''=cost% y'(0)=-1, y(0)=0% 轉(zhuǎn)換為一階微分方程 h=0.01;ht=10;t=0:h:ht;% y'=z1, z1'=y''=costy=zeros(1,length(t)

6、;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)

7、+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，求解區(qū)間0 25。MATLAB

8、程序：% 龍格庫塔歐拉方法解二階微分方程% 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

9、=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('