实验报告书8-常微分方程数值解.docVIP

东南大学《数学实验》报告 实验内容：常微分方程数值解 一 实验目的 自己编写常微分方程初值问题的常用算法，包括折线法、改进欧拉法、4阶龙格-库塔法（不允许直接使用ode45）,并用于对ODE模型的研究。 二 预备知识 (1)熟悉各种常用ODE数值算法原理 (2)了解各种算法的精度，熟悉ode45的用法 三 实验内容与要求 1.分别编写欧拉折线法、改进欧拉法和4阶龙格-库塔法通用算法 命令 欧拉法： function [t,x1,x2]=ForwardEuler(a,b,c,d,n) h=(b-a)/n; t=a:h:b; x1=[c zeros(1,n)]; x2=[d zeros(1,n)]; y=zeros(2,1); for i=1:n y=ODE(x1(i),x2(i)); x1(i+1)=x1(i)+h*y(1); x2(i+1)=x2(i)+h*y(2); end end 改进欧拉法： function [t,x1,x2]=ModifiedEuler(a,b,c,d,n) h=(b-a)/n; t=a:h:b; x1=[c zeros(1,n)]; x2=[d zeros(1,n)]; for i=1:n y=ODE(x1(i),x2(i)); yn1=x1(i)+h*y(1); yn2=x2(i)+h*y(2); dx=ODE(yn1,yn2); x1(i+1)=x1(i)+(h/2)*(y(1)+dx(1)); x2(i+1)=x2(i)+(h/2)*(y(2)+dx(2)); end end 4阶龙格-库塔法： function [t,x1,x2]=LK4(a,b,c,d,n) h=(b-a)/n; t=a:h:b; x1=[c zeros(1,n)]; x2=[d zeros(1,n)]; for i=1:n-1 k1=ODE(x1(i),x2(i)); xk2=ODE(x1(i)+h/2,x2(i)+h/2*k1(1)); yk2=ODE(x1(i)+h/2,x2(i)+h/2*k1(2)); k2=[xk2(1) yk2(2)]; xk3=ODE(x1(i)+h/2,x2(i)+h/2*k2(1)); yk3=ODE(x1(i)+h/2,x2(i)+h/2*k2(2)); k3=[xk3(1) yk3(2)]; xk4=ODE(x1(i)+h,x2(i)+h*k3(1)); yk4=ODE(x1(i)+h,x2(i)+h*k3(2)); k4=[xk4(1) yk4(2)]; x1(i+1)=x1(i)+h/6*(k1(1)+2*k2(1)+2*k3(1)+k4(1)); x2(i+1)=x2(i)+h/6*(k1(2)+2*k2(2)+2*k3(2)+k4(2)); end （2）用上述三种算法求解Lotka-Volterra模型（参数自行确定），并比较各种算法的计算精度 命令 结果 Lotka-Volterra模型: function dx=ODE(x1,x2) dx=zeros(2,1); dx(1)=x1*(1-0.1*x2); dx(2)=x2*(-0.5+0.02*x1); 主程序： [t,x1,x2]=ForwardEuler(0,15,25,2,150); plot(t,x1,-,t,x2,*) [t,x1,x2]=ModifiedEuler(0,15,25,2,150); plot(t,x1,-,t,x2,*) [t,x1,x2]=LK4(0,15,25,2,150); plot(t,x1,-,t,x2,*) function dx=shier1(t,x) dx=zeros(2,1); dx(1)=x(1)*(1-0.1*x(2)); dx(2)=x(2)*(-0.5+0.02*x(1)); [t,x]=ode45(shier1,[0 15],[25 2]); plot(t,x(:,1),-,t,x(:,2),*) 向前欧拉法： 改进欧拉法： 4阶龙格-库塔法： Ode45求得标准： 改进欧拉法精度最高，其他两种偏差较大。 （2）从下面两个问题中选作一个： 2－1：建立酒后驾驶模型教材(p131-132)，针对题目中的数据，分别用自己编写的龙格库塔法和Matlab内置的ode45命令进行研究。 2－2：自行收集足够的中国人口数