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

1、东南大学数学实验报告实验内容：常微分方程数值解一 实验目的自己编写常微分方程初值问题的常用算法，包括折线法、改进欧拉法、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

2、:n y=ODE(x1(i),x2(i); x1(i+1)=x1(i)+h*y(1); x2(i+1)=x2(i)+h*y(2);endend改进欧拉法：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/

3、2)*(y(2)+dx(2);endend4阶龙格-库塔法：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)

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

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

6、s(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）从下面两个问题中选作一个：21：建立酒后驾驶模型教材(p131-132)，针对题目中的数据，分别用自己编写的龙格库塔法和Matlab内置的ode45命令进行研究。22：自行收集足够的中国人口数据（可以是全国

7、或某个省市的数据），根据这些数据确定模型中可能用到的参数，并分别用指数增长模型和阻滞增长对人口情况做预报，并将理论和实际数据的进行比教。20002001200220032004200520062007200820094140 4186 4222 4254 4284 4311 4339 4368 4400 4432程序：clearclcdigits(2);year=2000:1:2009;RealPopulation=4140 4186 4222 4254 4284 4311 4339 4368 4400 4432;x0=4140;x1=4186;r=log(x1/x0);ExpPopu=x0.

8、*exp(r.*(year-year(1);errorEXP=100.*abs(ExpPopu-RealPopulation)./RealPopulation;rl=0.03402;xm=13040;LogPopu=xm./(1+(xm/x0-1).*exp(-r.*(year-year(1);errorLog=100.*abs(LogPopu-RealPopulation)./RealPopulation;plot(year,RealPopulation,'b',year,ExpPopu,'r');title('江西人口（Malthus）');figure;plot(year,100,'r*',year,errorEXP,'b');title('Malthus模型预测江西人口相对误差');figure;plot(year,RealPopulation,'b',year,Log