2013春数学实验基础 实验报告(1)常微分方程

年级、专业

姓名 学号 名单序号 实验时间 2013年 3月 日 使用设备、软件 PC, MATLAB 注： 实验报告的最后一部分是实验小结与收获

实验一 常微分方程

8 / 9 2013春 数学实验 实验一 常微分方程 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 法','解析解');

年级、专业 姓名 学号 名单序号

实验时间 2013年 3月 日 使用设备、软件 PC, MATLAB

注： 实验报告的最后一部分是实验小结与收获

实验一 常微分方程

8 / 9

2013春 数学实验 实验一 常微分方程 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 处的切线斜率等于22,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');

年级、专业 姓名 学号 名单序号 实验时间 2013年 3月 日 使用设备、软件 PC, MATLAB 注： 实验报告的最后一部分是实验小结与收获

实验一 常微分方程

8 / 9

2013春 数学实验 实验一 常微分方程

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,'*');

年级、专业 姓名 学号 名单序号

实验时间 2013年 3月 日 使用设备、软件 PC, MATLAB

注： 实验报告的最后一部分是实验小结与收获

实验一 常微分方程

8 / 9

2013春 数学实验 实验一 常微分方程

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 '，由此建立微分方程为

年级、专业 姓名 学号 名单序号

实验时间 2013年 3月 日 使用设备、软件 PC, MATLAB

注： 实验报告的最后一部分是实验小结与收获

实验一 常微分方程

8 / 9

2013春 数学实验 实验一 常微分方程 ?????=-?='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 ;

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

(2) 厂家问：要做多少时间广告，可使市场购买率达到80%？

t_min=min(find(x>0.8));

t(t_min)

ans =

8.8543

结果：大约8.8543个单位的时间后，可使市场购买率达到80%。

6. (肿瘤生长) 肿瘤大小V 生长的速率与V 的a 次方成正比，其中a 为形状参数，10≤≤a ;而其比例系数K 随时间减小，减小速率又与当时的K 值成正比，比例系数为环境参数b 。设某肿瘤参数a=1, b=0.1, K 的初始值为2，V 的初始值为1。问

(1)此肿瘤生长不会超过多大？

年级、专业 姓名 学号 名单序号 实验时间 2013年 3月 日 使用设备、软件 PC, MATLAB 注： 实验报告的最后一部分是实验小结与收获

实验一 常微分方程

8 / 9 2013春 数学实验 实验一 常微分方程 由已知条件建立微分方程：

???????==?-='≤≤?='2

…… 此处隐藏：3911字，全部文档内容请下载后查看。喜欢就下载吧 ……