数学实验之微分方程

数学实验之微分方程第一页，共四十七页，编辑于2023年，星期六

MATLAB软件的实现引例：增长与传播模型

实验内容

微分方程模型与方法

解析解

数值解主要内容第二页，共四十七页，编辑于2023年，星期六

1)掌握微分方程求解的三种解法：解析法、数值解法以及图形表示解的方法;2)掌握求微分方程数值解的欧拉方法,了解龙格---库塔方法的思想;3）学会使用MATLAB软件求解析解、数值解和图形解；

4）通过范例学习怎样建立微分方程模型和分析问题的思想；实验目的返回第三页，共四十七页，编辑于2023年，星期六首先研究某些物种增长的一阶微分方程.其中x(t)表示一种给定物种在时刻t的总数;r(t,x)表示该物种的净增长率;进一步假定r(t,x)是常数,于是得到简化模型：解：引例:群体增长模型第四页，共四十七页，编辑于2023年，星期六t=0:12;x=2*exp(0.4*t);plot(t,x,'r:o')思考:该模型是否符合物种的自然增长规律?第五页，共四十七页，编辑于2023年，星期六

该模型用于预测地球上的人口总数,例如1961年地球上的人口总数为3,060,000,000，又假定人口总数以2%的速度增长(净增长率)，故检验：1700~1961年期间，实际数字与理论数值吻合！预测：p(2510)=2,000,000(亿)；

p(2670)=36,000,000(亿)；星球的总表面积约18,600,000亿ft2，80%被水覆盖，到2510年，每人占地面积只有9.3

ft2,到2670年，每人占地面积只有0.5167ft2。第六页，共四十七页，编辑于2023年，星期六分析以上模型,并修改模型.一般地，b<<r，若x不是很大，则可以忽略-bx2。相反,当x很大时,-bx2不能忽略,它将减慢群体迅速的增长率.解出结果：竞争项第七页，共四十七页，编辑于2023年，星期六分析解

称该曲线为逻辑曲线或S形曲线。符合很多实际情况。

物种发展到一定时期,会受到其它因素的制约,增长速度渐趋缓慢.第八页，共四十七页，编辑于2023年，星期六引例：新技术的传播模型

一项新技术如何在同类企业中推广?传播速度怎样？哪些因素影响传播速度？第九页，共四十七页，编辑于2023年，星期六1）记p(t)表示t时刻采用了该技术的企业数，并假定它连续可微；2）该项技术具有成效,值得传播；模型假设3）当t=0时，有一个企业已经采用该技术，共有N个同类企业数；4）在[t,t+△t]内，△p(t)与p(t)(N-p(t))△t成正比；第十页，共四十七页，编辑于2023年，星期六由条件4），得到模型建立新技术传播模型

某行业内采纳一项新技术的企业数必然依赖于该技术的效益L和所需要的投资资金s，一般情况下,一个企业的日传播率r=f(L,s).但这里假定r为常数.第十一页，共四十七页，编辑于2023年，星期六模型求解计算结果：逻辑曲线①P(t)>0，并且t→+∞时，p(t)→N;②对任意的t，p(t)<N;③当p=N/2时，dp/dt达到最大；结果分析第十二页，共四十七页，编辑于2023年，星期六模型检验

该模型需要进行实证分析，检验是否符合实际情况，检验模型假设是否合理。

例如,考察过去的内燃机火车头、集中化的交通控制、车厢减速器等三种技术从一个企业推广到另一些企业的速率问题，如下图所示。第十三页，共四十七页，编辑于2023年，星期六理论曲线N=25t0=1925实测数据p0=1r1=4%r2=3.1%r3=1.5%第十四页，共四十七页，编辑于2023年，星期六模型修正

随着通讯能力的提高和大众媒介的普及，广告将对这项新技术的传播起到一定的作用。因此，在上述模型中增加一个因素——广告效应的作用。即修改模型如下：计算结果：仍然是逻辑曲线，与实际情况应该是更加吻合。第十五页，共四十七页，编辑于2023年，星期六1、上述模型哪些假设需要修改？

2、流行性感冒的传播模型；思考:第十六页，共四十七页，编辑于2023年，星期六定义：含有导数的方程称为微分方程。如上例，

f(x,V(x),V’(x))=0微分方程模型1、微分方程的一般形式：

F(x,y,y’,…,y(n))=0

隐式或

y(n)=f(x,y,y’,…,y(n-1))

显式n阶第十七页，共四十七页，编辑于2023年，星期六特殊情形：2、一阶微分方程组的一般形式：1阶初始条件：y(x0)=y0微分方程模型第十八页，共四十七页，编辑于2023年，星期六满足一定的初始条件，称为偏微分方程。微分方程模型第十九页，共四十七页，编辑于2023年，星期六微分方程求解方法简介③图形解xyo①简单的微分方程②复杂、大型的微分方程返回①

解析解

y=f(x)②

数值解

(xi,yi)欧拉方法改进欧拉方法

梯形法龙格-库塔法第二十页，共四十七页，编辑于2023年，星期六MATLAB软件实现解析解dsolve('eqn1','eqn2',…,'c1',…,'var')微分方程组初值条件自变量名注意：①y'Dy，y''D2y②自变量名可以省略，默认变量名‘t’。第二十一页，共四十七页，编辑于2023年，星期六例①输入：y=dsolve('Dy=1+y^2')y1=dsolve('Dy=1+y^2','y(0)=1','x')输出：y=tan(t-C1)（通解，一簇曲线）

y1=tan(x+1/4*pi)（特解，一条曲线）MATLAB软件实现第二十二页，共四十七页，编辑于2023年，星期六例②常系数的二阶微分方程y=dsolve('D2y-2*Dy-3*y=0','x')y=dsolve('D2y-2*Dy-3*y=0','y(0)=1,Dy(0)=0','x')输入：x=dsolve('D2x-(1-x^2)*Dx+x=0','x(0)=3,Dx(0)=0')上述两例的计算结果怎样？由此得出什么结论？例③非常系数的二阶微分方程例③无解析表达式！第二十三页，共四十七页，编辑于2023年，星期六x=dsolve('(Dx)^2+x^2=1','x(0)=0')例④非线性微分方程x=[sin(t)][-sin(t)]若欲求解的某个数值解，如何求解？t=pi/2;eval(x)MATLAB软件实现第二十四页，共四十七页，编辑于2023年，星期六输入：[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y')[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y','x(0)=0,y(0)=1')例④输出：（li3.m）MATLAB软件实现返回第二十五页，共四十七页，编辑于2023年，星期六数值解1、欧拉法2、龙格—库塔法数值求解思想：(变量离散化)

引入自变量点列{xn}→{yn}，在x0x1x2…xn…上求y(xn)的近似值yn.通常取等步长h，即xn=x0+n×h，或

xn

=xn-1+h，（n=1,2,…)。研究常微分方程的数值解法是十分必要的。解:y=y(x)第二十六页，共四十七页，编辑于2023年，星期六1)向前欧拉公式:（y’=f(x,y)）

y(xn+1)

y(xn)+hf(xn,y(xn))(迭代式)

yn+1

yn+hf(xn,yn)(近似式)

特点：f(x,y)取值于区间[xn,xn+1]的左端点.1、欧拉方法在小区间[xn,

xn+1]上用差商代替微商(近似),第二十七页，共四十七页，编辑于2023年，星期六

yn+1

yn

+hf(xn+1,yn

+1)特点：①f(x,y)取值于区间[xn,xn+1]的右端点.②非线性方程,称‘隐式公式’。yn+1

=

yn+hf(xn,yn)2)向后欧拉公式方法：迭代（y’=f(x,y)）x=[];y=[];x(1)=x0;y(1)=y0;forn=1:kx(n+1)=x(n)+n*h;y(n+1)=

y(n)+h*f(x(n),y(n));(向前）end,x,y,1、欧拉方法第二十八页，共四十七页，编辑于2023年，星期六例1

观察向前欧拉、向后欧拉算法计算情况。与精确解进行比较。误差有多大？解：1）解析解：y=x+e-xy=dsolve('Dy=-y+x+1','y(0)=1','x')1、欧拉方法第二十九页，共四十七页，编辑于2023年，星期六2）向前欧拉法：

yn+1=yn+h(-yn+xn+1)=(1-h)yn+hxn+h3）向后欧拉法：

yn+1=yn+h(-yn+1+xn+1+1)

转化yn+1=(yn+hxn+1+h)/(1+h)y’=f(x,y)=-y+x+1;1、欧拉方法第三十页，共四十七页，编辑于2023年，星期六x1(1)=0;y1(1)=1;y2(1)=1;h=0.1;%(died.m)fork=1:10x1(k+1)=x1(k)+h;y1(k+1)=(1-h)*y1(k)+h*x1(k)+h;y2(k+1)=(y2(k)+h*x1(k+1)+h)/(1+h);endx1,y1,y2,%（y1——向前欧拉解，y2——向后欧拉解）x=0:0.1:1;y=x+exp(-x)%（解析解）plot(x,y,x1,y1,'k:',x1,y2,'r--')1、欧拉方法第三十一页，共四十七页，编辑于2023年，星期六x精确解向前欧拉向后欧拉01110.11.004811.00910.21.01871.011.02640.31.04081.0291.05130.41.07031.05611.08300.51.10651.09051.12090.61.14881.13141.16450.71.19661.17831.21320.81.24931.23051.26650.91.30661.28741.324111.36791.34871.3855（1）步长h=0.1的数值解比较表计算结果第三十二页，共四十七页，编辑于2023年，星期六（2）步长h=0.01的数值解比较表x精确解向前欧拉向后欧拉01110.11.00481.00441.00530.21.01871.01791.01950.31.04081.03971.04190.41.07031.06901.07170.51.10651.10501.10800.61.14881.14721.15040.71.19661.19481.19830.81.24931.24751.25110.91.30661.30471.308411.36791.36601.3697结论：显然迭代步长h的选取对精度有影响。第三十三页，共四十七页，编辑于2023年，星期六图形显示

有什么方法可以使精度提高？向后欧拉法返回第三十四页，共四十七页，编辑于2023年，星期六

对方程y’=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：2、使用数值积分即梯形法：第三十五页，共四十七页，编辑于2023年，星期六梯形公式

改进欧拉公式yn+hf(xn,yn)返回第三十六页，共四十七页，编辑于2023年，星期六

以例1为例，用改进欧拉公式编程计算，再与精确解的比较。yn+1=yn+(h/2)*[(-yn+xn+1)+(-yn+1+xn+1+1)]=yn+(h/2)*[(-yn+xn+1)-(yn+h*(-yn+xn+1))+xn+1+1]=yn+(h/2)*[(1-h)*xn+xn+1+2-h+(h-2)*yn]died1.m使用改进欧拉公式的例第三十七页，共四十七页，编辑于2023年，星期六x精确解向前欧拉向后欧拉改进欧拉011110.11.004811.00911.0050.21.01871.011.02641.0190.31.04081.0291.05131.04120.41.07031.05611.08301.07080.51.10651.09051.12091.10710.61.14881.13141.16451.14940.71.19661.17831.21321.19720.81.24931.23051.26651.25000.91.30661.28741.32411.307211.36791.34871.38551.3685步长

h=0.1的数值解比较表使用改进欧拉公式的例第三十八页，共四十七页，编辑于2023年，星期六3、使用泰勒公式

以此方法为基础，有龙格-库塔法、线性多步法等方法。4、数值公式的精度

当一个数值公式的截断误差可表示为O（hk+1）时（k为正整数，h为步长），称它是一个k阶公式。

k越大，则数值公式的精度越高。

欧拉法是一阶公式，改进的欧拉法是二阶公式。龙格-库塔法有二阶公式和四阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。返回第三十九页，共四十七页，编辑于2023年，星期六[t，x]=solver（’f’,ts,x0）ode23

ode45由待解方程写成的m-函数文件ts=[t0，tf]，t0、tf为自变量的初值和终值函数的初值ode23：组合的2/3阶龙格-库塔算法ode45：运用组合的4/5阶龙格-库塔算法自变量值函数值Matlab软件计算数值解第四十页，共四十七页，编辑于2023年，星期六1）首先建立M-文件（weif.m）

functionf=weif(x,y)f=-y+x+1;2）求解：[x，y]=ode23(‘weif’,[0,1],1)3)作图形:plot(x,y,‘r’);4)与精确解进行比较

holdonezplot(‘x+exp(-x)’,[0,1])例1

y’=-y+x+1，y(0)=1标准形式：y’=f(x,y)范例第四十一页，共四十七页，编辑于2023年，星期六1、在解n个未知函数的方程组时，x0和x均为n维向量，m-函数文件中的待解方程组应以x的分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意:第四十二页，共四十七页，编辑于2023年，星期六注意1：1、建立M文件函数

functionxdot=fun(t,x,y)xdot=[f1(t,x(t),y(t))；

f2(t,x(t),y(t))]；2、数值计算（执行以下命令）

[t,x,y]=ode23(‘fun',[t0,tf],[x0,y0])注意：执行命令不能写在M函数文件中。xd(1)=f1(t,x(t),y(t));xd(2)=f2(t,x(t),y(t));xdot=xd’;%