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

)0(1)0()

()(1

0,)()()(K V t K b t K a t V t K t V a

%先编写上述函数的ex6_1fun.m 文件

function f=ex6_1fun(t,y)

f(1)=y(2).*y(1);

f(2)=-0.1*y(2);

f=f(:);

%画出其图像，并求最大的肿瘤大小V

[t1,y1]=ode45(@ex6_1fun,[0,100],[1,2]);

plot(t1,y1(:,1),'r*');

max(y1(:,1))

8

ans =

4.8567e+008

%故肿瘤生长不会超过8108567.4?

(2)过多长时间肿瘤大小翻一倍？

%从图像可以看出，大约在（0，1）内，肿瘤大小翻一倍，以此求解

[t2,y2]=ode45(@ex6_1fun,[0 1],[1;2]);

t_min=min(find(y2(:,1)>2));

t2(t_min)

ans =

0.3750

答：大约0.4个单位时间后肿瘤大小翻一倍。

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

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

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

实验一常微分方程

(3)何时肿瘤生长速率由递增转为递减？

dv=y1(:,2).*y1(:,1);

[Vn,tn]=max(dv);

t1(tn)

ans =

29.5466

答：大约30个单位时间后肿瘤生长速率由递增转为递减。

(4)若参数a=2/3呢？

%重新编写微分方程ex6_4fun，并依次计算上述三个问题

function f=ex6_4fun(t,y)

f(1)=y(2)*y(1).^(2/3);

f(2)=-0.1*y(2);

f=f(:);

%画图像，肿瘤生长不会超过450.7959

%求多长时间肿瘤大小翻一倍，大约为0.4

ans =

0.4000

%何时肿瘤生长速率由递增转为递减，大约为10

ans =

9.5718

选做题：

1.(生态系统的振荡现象)第一次世界大战中，因为战争很少捕鱼，按理战后应能捕到更多

的鱼才是。可是大战后，在地中海却捕不到鲨鱼，因而渔民大惑不解。

令x1为鱼饵的数量，x2为鲨鱼的数量，t为时间。微分方程为

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

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

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

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

实验一 常微分方程

8 / 9

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

式中a 1, a 2, b 1, b 2都是正常数。第一式鱼饵x 1的增长速度大体上与x 1成正比，即按a 1x 1比率增加, 而被鲨鱼吃掉的部分按b 1x 1x 2的比率减少；第二式中鲨鱼的增长速度由于生存竞争的自然死亡或互相咬食按a 2x 2的比率减少，但又根据鱼饵的量的变化按b 2x 1x 2的比率增加。对a 1=3, b 1=2, a 2=2.5, b 2=1, x 1(0)=x 2(0)=1求解。画出解曲线图和相轨线图，可以观察到鱼饵和鲨鱼数量的周期振荡现象。

%首次编写上述微分方程的ex_7fun.m 函数

function f=ex_7fun(t,x)

a1=3;b1=2;a2=2.5;b2=1;

f(1)=x(1)*(a1-b1*x(2));

f(2)=-x(2)*(a2-b2*x(1));

f=f(:);

%画出解曲线图和相轨线图

[t,x]=ode45(@ex_7fun,[0,4],[1;1]);

subplot(1,2,1);

plot(t,x(:,1),'r',t,x(:,2),'k:');

legend('x1','x2');

title('解曲线');

subplot(1,2,2);

plot(x(:,1),x(:,2));

title('相轨线'); 0.51

1.5

2

2.5

3

3.5

4

4.5

55.5

x1x2

2.编写四阶Runge-Kutta 法程序

实验小结与收获：完成第二份实验报告后，我学习了Matlab 有关求解常微分方程（组）的指令，了解了什么是Euler 法和刚性方程组以及如何求解，并了解了ode45、ode15s 、Euler 等函数间的异同；某些无法求得解析解的方程可以通过计算数值解了解它们的性质，有助于我们的学习；还学习了通过建立数学模型解决实际问题，将理论应用于实际，锻炼了我们逻辑思维能力和抽象概括能力。

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

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

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

实验一常微分方程

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

本文来源：https://www.bwwdw.com/article/2hve.html