實(shí)驗(yàn)報(bào)告書8-常微分方程數(shù)值解

1、東南大學(xué)數(shù)學(xué)實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)內(nèi)容：常微分方程數(shù)值解一 實(shí)驗(yàn)?zāi)康淖约壕帉懗Ｎ⒎址匠坛踔祮栴}的常用算法，包括折線法、改進(jìn)歐拉法、4階龍格-庫塔法（不允許直接使用ode45）,并用于對ODE模型的研究。 二 預(yù)備知識 (1)熟悉各種常用ODE數(shù)值算法原理(2)了解各種算法的精度，熟悉ode45的用法三 實(shí)驗(yàn)內(nèi)容與要求1.分別編寫歐拉折線法、改進(jìn)歐拉法和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改進(jìn)歐拉法：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模型（參數(shù)自行確定），并比較各種算法的計(jì)算精度命令結(jié)果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),'*')向前歐拉法：改進(jìn)歐拉法：4階龍格-庫塔法：Ode45求得標(biāo)準(zhǔn)：改進(jìn)歐拉法精度最高，其他兩種偏差較大。（2）從下面兩個(gè)問題中選作一個(gè)：21：建立酒后駕駛模型教材(p131-132)，針對題目中的數(shù)據(jù)，分別用自己編寫的龍格庫塔法和Matlab內(nèi)置的ode45命令進(jìn)行研究。22：自行收集足夠的中國人口數(shù)據(jù)（可以是全國

7、或某個(gè)省市的數(shù)據(jù)），根據(jù)這些數(shù)據(jù)確定模型中可能用到的參數(shù)，并分別用指數(shù)增長模型和阻滯增長對人口情況做預(yù)報(bào)，并將理論和實(shí)際數(shù)據(jù)的進(jìn)行比教。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模型預(yù)測江西人口相對誤差');figure;plot(year,RealPopulation,'b',year,Log