龙格库塔解一阶微分方程c语言,四阶龙格库塔法解微分方程

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

精品资料 欢迎下载 四阶龙格库塔法解微分方程 一、四阶龙格库塔法解一阶微分方程 ①一阶微分方程： ，初始值y(0)=0，求解区间[0 10]。 MATLAB程序： %%%%%%%%%%% 四阶龙哥库塔法解一阶微分方程 %%%%%%%%%%% y =cost %%%%%%%%%%% y(0)=0, 0≤t≤10，h=0.01 %%%%%%%%%%% 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 ); 图1 ②一阶微分方程： ，初始值x(1)=2，求解区间[1 3]。 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; 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 =cost y=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_