春数学实验基础实验报告常微分方程

1. 分别用Euler 法和ode45解下列常微分方程并与解析解比较： (1)30,1)0(,<<=+='x y y x y

编写Euler 法的matlab 函数，命名为euler.m

function [t,y]=euler(odefun,tspan,y0,h) t=tspan(1):h:tspan(2); y(1)=y0;

for i=1:length(t)-1

y(i+1)=y(i)+h*feval(odefun,t(i),y(i)); end t=t'; y=y';

下面比较三者的差别： % ode45

odefun=inline('x+y','x','y'); [x1,y1]=ode45(odefun,[0,3],1); plot(x1,y1,'ko');pause hold on ; % Euler·¨

[x2,y2]=euler(odefun,[0,3],1,0.05); plot(x2,y2,'r+');pause hold on ;

% 解析解

y0=dsolve('Dy=t+y','y(0)=1'); ezplot(y0,[0,3]);pause hold off ;

legend('ode45','euler 法','解析解');

Euler 法只有一阶精度，所以实际应用效率比较差，而ode45的效果比较好，很接近真实值。 (2) 2

0.01()2sin(),(0)0,(0)1,05y y y t y y t ''''-+===<<

先写M 文件ex1_2fun.m

function f=ex1_2fun(t,y) f(1)=y(2);

f(2)=0.01*y(2).^2-2*y(1)+sin(t); f=f(:); % ode45

[t1,y1]=ode45(@ex1_2fun,[0,5],[0;1]); plot(t1,y1(:,1),'ko');

% 解析解

s=dsolve('D2y-0.01*(Dy)^2+2*y=sin(t)','y(0)=0','Dy(0)=1','t')

s =

[ empty sym ]

%由此可知该微分方程无解析解

2. 求一通过原点的曲线，它在(,)x y 处的切线斜率等于2

2,0 1.57.x y x +<<若x 上限增为1.58，1.60会发生什么？

odefun=inline('2*x+y^2','x','y'); subplot(1,4,1);

[x1,y1]=ode45(odefun,[0,1.57],0); plot(x1,y1,'r*'); title('上限1.57'); subplot(1,4,2);

[x2,y2]=ode45(odefun,[0,1.58],0); plot(x2,y2,'bo'); title('上限1.58'); subplot(1,4,3);

[x3,y3]=ode45(odefun,[0,1.6],0); plot(x3,y3,'k'); title('上限1.60'); subplot(1,4,4);

plot(x1,y1,'r*');hold on ; plot(x2,y2,'bo');hold on ; plot(x3,y3,'k');hold off ;

legend('上限1.57','上限1.58','上限1.60');

13

14

14

结论：随着x 上界的增加，解趋于无穷大。 3. 求解刚性方程组：

500,1)0(,5.025.100075.999,1)0(,5.075.99925.10002

2121211

<<⎩⎨

⎧-=+-='=++-='x y y y y y y y y 先写M 函数ex3fun.m

function f=ex3fun(t,y)

f(1)=-1000.25*y(1)+999.75*y(2)+0.5; f(2)=999.75*y(1)-1000.25*y(2)+0.5; f=f(:);

%作图

[t,y]=ode15s(@ex3fun,[0,50],[1,-1]); plot(t,y,'*');

4. (温度过程)夏天把开有空调的室内一支读数为20℃的温度计放到户外，10分钟后读2

5.2℃, 再过10分钟后读数28.32℃。建立一个较合理的模型来推算户外温度。 设：t 时刻温度计的读数为T ，户外温度为c ，T 的增速与室内外温差（c-T ）成正比，由此建立微分方程，其中k 为比例系数

⎩⎨

⎧=-='20

)0()

(T T c k T %首先，计算解析解

>> y=dsolve('DT=k*(c-T)','T(0)=20','t')

y =

c - (c - 20)/exp(k*t)

%又已知32.28)20(,2.25)10(==T T ，用非线性最小二乘拟合该函数，调用lsqcurvefit 命令：

fun=inline('c(1)-(c(1)-20)./exp(c(2)*t)','c','t'); lsqcurvefit(fun,[30 1],[10 20],[25.2 28.32])

ans =

33.0000 0.0511 即户外温度c=33，比例系数k=0.0511

5. (广告效应)某公司生产一种耐用消费品，市场占有率为5%时开始做广告，一段时间的市场跟踪调查后，该公司发现：单位时间内购买人口百分比的相对增长率与当时还没有买的百分比成正比，且估得此比例系数为0.5。

(1) 建立该问题的数学模型，并求其数值解与模拟结果作以比较； 设：市场占有率为x ，市场占有率的增长速度为x '，则相对增长率为

x

x '

，由此建立微分方程为 ⎪⎩⎪⎨⎧=-⨯='

05

.0)0()

1(5.0x x x

x %首先，计算解析解

>> s=dsolve('Dx=(0.5*(1-x))*x','x(0)=0.05','t') s=

1/(exp(log(19) - t/2) + 1)

%再调用ode45计算数值解，并作图比较解析解与数值解的区别： odefun=inline('(0.5*(1-x))*x','t','x'); [t,x]=ode45(odefun,[0,20],0.05); plot(t,x,'r*');hold on ; ezplot(s,[0 20]);hold off ;