最新四阶龙格库塔法解微分方程.docxVIP

精品 精品 - - 可编辑 - 精品 精品 - - 可编辑 - 四阶龙格库塔法解微分方程 一、四阶龙格库塔法解一阶微分方程 ①一阶微分方程： y cost ，初始值 y(0)=0 ，求解区间 [0 10] MATLAB 程序： %%%%%%%%%%% 四阶龙哥库塔法解一阶微分方程 %%%%%%%%%%% y=cost %%%%%%%%%%% y(o)=o, 0 t10 , h=o.oi %%%%%%%%%%% y=sint h=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 ); 精品 精品 - -可编辑- 精品 精品 - - 可编辑 - Approximation Exact 图1 2 2 ②一阶微分方程：X t*) x /t ,初始值x(1)=2 ,求解区间[1 3] MATLAB 程序： %%%%%%%%%%% 四阶龙哥库塔法解微分方程 %%%%%%%%%%% x(t)=(t*x-xA2)/tA2 %%%%%%%%%%% x(1)=2, 1 t 3, h=1/128 %%%%%%%%%%% 精确解：x(t)=t/(0.5+l nt) %%%%% 步长 %%%%% 步长 tf=3; t=1:h:tf; x=zeros(1,length(t)); x(1)=2; %%%%% 初始值 F_tx=@(t,x)(t.*x-x.A2)./t.A2; 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); 二、四阶龙格库塔法解二阶微分方程 ①二阶微分方程：y cost，初始值y(0)=0 ,y(0)=-1，求解区间[0 10] MATLAB 程序： %%%%%%%%% 龙格库塔欧拉方法解二阶微分方程 %%%%%%%%% y=cost vh(l)z oh(lm a)llo)ue「L)so」① ZHZ a)llo)ue「L)so」① ZUA ls8Hb-N-N女 ％%%%%%%%% AflloL OLILU =0.0 丄 魁40箕——蔓牌 ％%%%%%%%% OOMO)-A %%%%%%%%% 龔 -Ik- ◎_J*LI + ( 一Ncox*ll+(一 M-ll+§)6寸_1 心_1*1|+(一Ncox*ll+(一 M-ll+slr寸 x c?_J*2a+(DZ7x*Qll+(_M7a+§)6Hco_J c?_J*2a+(DZ7x*Qll+(一 M0ll+§1rcox 二—lda+aNLX*QLI+(_M7a+§)6H2_J 二_J*0a+(_NLX*QLI+(_M0LI+§1r2x gzeA-slrLX ?e￡6u①三l上」0」 a)u_s)(NA」)?6 _(Z)(NA」)◎上 - -可编辑- - -可编辑- 精品 y(i+1)=y(i)+h/6*(K1+2*K2+2*K3+K