第八章微分方程

数学实验第八章微分方程

8.1引例解:设所求曲线方程为y=y(x),则又y(x)应满足条件x=1时y=2，因此可解得C=1两端对x积分，得例1求xOy平面上过点(1,2)的一条曲线方程，使其在任意点处的斜率为2x。故所求曲线方程为y=x2+1。例2列车上以20米/秒的速度行驶，当制动时，列车获得加速度-0.4米/秒2，问开始制动后多长时间能停住，以及列车在这段时间行驶了多少路程？等式两边同时积分，得解:设列车开始制动后t秒内行驶了s米，则再积分得利用条件t=0,s=0,ds/dt=20，可以解得C1=20,C2=0,因此利用条件列车停住时速度为零，由解得:s=500(米)代入及C1=20可得:t=50(秒)微分方程:例如8.2微分方程的基本概念含有未知函数的导数或微分的方程常微分方程(ODE)(OrdinaryDifferentialEquations)仅含一个自变量的微分方程.微分方程的阶:

微分方程中出现的未知函数的最高阶导数的阶数.线性与非线性微分方程:方程中关于未知函数及其各阶导数均是一次的,则称为线性微分方程.线性微分方程非线性微分方程微分方程的解:代入微分方程能使方程成为恒等式的函数称为微分方程的解。微分方程解的分类：(1)通解:微分方程的解中含有任意常数,且任意常数的个数与微分方程的阶数相同。(2)特解:确定了通解中任意常数以后的解。解的图象(积分曲线):

微分方程的解曲线。初始条件:

用来确定任意常数的条件。考虑引例1，对于微分方程函数y=x2+C是其通解，y=x2+1是其满足初始条件x0=1，y(x0)=2的特解。微分方程的(部分)积分曲线(0≤C≤2)见右图，其中红色曲线为满足初始条件的特解曲线。考虑引例2，对于微分方程函数s=-0.2t2+C1t+C2是其通解,s=-0.2t2+20t是其满足条件t0=0，s(t0)=0,ds/dt(t0)=20的特解.满足初始条件的积分(特解)曲线见下图。路程曲线速度曲线

微分方程是研究函数变化规律的有力工具，在科技、工程技术、经济管理以及生态、环境、人口、交通等各个领域中有着十分广泛的应用。

很多实际问题的数学模型都可以用微分方程来表示，但建立微分方程只是解决问题的第一步，通常在此基础上还需要求出微分方程的解并对实际问题进行解释和加以检验。

求解微分方程一般可分为求微分方程的解析解和数值解。微分方程的解

解析解可以用一个确切的代数表达式来表示的解。显然微分方程的解析解对微分方程的分析和应用都是很方便的。但是存在解析解的微分方程仅仅局限于一些特殊类型，而绝大多数微分方程都是求不出解析解的。

数值解求解不存在解析解的微分方程时必须借助数值解法去求出微分方程的数值解(近似解)。因此，数值解法是求解微分方程的一个十分重要的方法。ODE的解很少能用初等函数及其不定积分的组合表示。例如：方程不能表示为初等函数，故得不到精确解的解是难以求积方程8.3常微分方程数值解法介绍(一)常微分方程数值解的定义在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解析解。而在实际上，对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式。因此，研究常微分方程的数值解法是十分必要的。对一阶微分方程

dy/dx=f(x,y),y(x0)=y0主要方法(1)欧拉方法(向前欧拉公式，向后欧拉公式，梯形公式，改进的欧拉公式)(2)龙格-库塔方法(二阶，三阶，四阶，五阶)在自变量x的取值点列{xn}，根据一定的原理(主要是用差商代替导数)，利用迭代方法求出y(xn)的近似值yn(二)

建立数值解法的一些途径1、用差商代替导数若步长h较小，则有故有公式：此即欧拉法。设xi+1-xi=h,(i=0,1,…,n-1),可以用以下离散化方法求解微分方程dy/dx=f(x,y),y(x0)=y02、使用数值积分对方程y'=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：故有公式：实际应用时，与欧拉公式结合使用：此即改进的欧拉法。3、使用泰勒公式以此方法为基础，有龙格-库塔法、线性多步法等方法。4、数值公式的精度当一个数值公式的截断误差y(xi)-yi可表示为O(hk+1)时(k为正整数，h为步长)，称它是一个k阶公式。k越大，则数值公式的精度越高。欧拉法是一阶公式，改进的欧拉法是二阶公式。龙格-库塔法有二至五阶公式。线性多步法有四阶外插公式和内插公式。8.4微分方程求解的MATLAB实现

一、求微分方程(组)的解析解求微分方程(组)解析解的MATLAB命令dsolve

1．求微分方程的通解格式：dsolve('eqn','v')

其中'eqn'

是微分方程表达式构成的字符串，'v'是微分方程表达式中自变量符号。'v'可以省略，省略时默认t为自变量。

在表达微分方程时,用字母D表示微分,D2,D3表示二,三阶微分。即在不指定自变量(默认t为自变量)时，Dy=dy/dt，D2y=d2y/dt2,…例8.1求微分方程dy/dx=1+y2的通解

dsolve('Dy=1+y^2')

ans=tan(t-C1)其中t为默认的自变量。如果希望以x为自变量，则需用命令

dsolve('Dy=1+y^2','x')

ans=tan(x-C1)注意：该命令中自变量的取法只要不与因变量相同即可。C1cos(x)+C2sin(x)---------------------1/2x例8.2求微分方程

x2y″+xy′+(x2-1/4)y=0的通解s=dsolve('x^2*D2y+x*Dy+(x^2-1/4)*y=0','x')

s=(C1*cos(x)+C2*sin(x))/x^(1/2)注意(1)若方程是右端等于0的形式，可以省略。(2)若不加自变量x，则将把x作为常数求解。(3)所得结果可以用命令pretty进行简化。例pretty(s)2.求微分方程的特解格式:dsolve('eqn','y(a)=b0,Dy(a)=b1,···','v')其中'eqn'是微分方程表达式构成的字符串，'y(a)=b0,Dy(a)=b1,···'

是微分方程的初始条件。'v'是微分方程表达式中自变量符号。例8.3求微分方程y''+4y'+29y=0满足条件y(0)=0,y’(0)=15的特解。dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')

ans=3*exp(-2*x)*sin(5*x)3.求解微分方程组

格式：dsolve('eqn1','eqn2',…)例8.4求下面微分方程组的通解。

dx/dt=3x+4y，dy/dt=-4x+3y

[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y')x=exp(3*t)*cos(4*t)*C1+exp(3*t)*sin(4*t)*C2y=-exp(3*t)*sin(4*t)*C1+exp(3*t)*cos(4*t)*C2注意：也可以用此命令求满足条件的特解

[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y','x(0)=1,y(0)=0')x=exp(3*t)*cos(4*t)y=-exp(3*t)*sin(4*t)二、求微分方程(组)的数值解注意求微分方程(组)数值解的MATLAB命令ode微分方程的数值解命令ode仅仅适用于一阶微分方程(组)。对高阶微分方程必须先将其化为一阶微分方程(组)再求解。高阶微分方程化为一阶微分方程组的方法

对任意高阶微分方程

x(n)=f(t,x,x',···,x(n-1))

令y1=xy2=y1'=x'y3=y2'=x''···yn=y'n-1=x(n-1)y1'=y2y2'=y3···yn-1'=ynyn'

=f(t,y1,y2,···,yn)

所以下面只讨论一阶微分方程组的数值解

一阶微分方程组的向量形式表示

令一阶微分方程组通常可用向量形式表示为dy/dt=f(t,y)

一阶微分方程组的向量形式表示

dy/dt=f(t,y)一阶线性微分方程组的向量形式表示

dy/dt=Ay其中求一阶微分方程组数值解的MATLAB命令ode

ode45,ode23,ode113,ode15s,ode23s五个命令的格式完全一致,常用的是前两个。格式：[T,Y]=ode45('fun',Ts,Y0)'fun'：由微分方程组写成的函数M文件的文件名；Ts=[t1,t2]：方程组解的取值区间；Y0

：方程组解的初始值(n维向量)输出T是m维列向量(取值从t1到t2)；Y是m×n维矩阵，其第i列是与T相应的yi的值。向量微分方程组dy/dt=f(t,y)的ode输出[T,Y]=ode45('fun',Ts,Y0,options)

例8.5求解线性微分方程组X'=AX，其中x1'

=x2x2'=4x1+3x2-4x3x3'=x1+2x2+x3

解：建立函数M文件exam1：functiondx=exam1(t,x)dx=zeros(3,1);dx(1)=x(2);dx(2)=4*x(1)+3*x(2)-4*x(3);dx(3)=x(1)+2*x(2)+x(3);取初始条件x(0)=[1;2;3],Ts=[0,20],利用ODE45命令求解如下[T,Y]=ode45('exam1',[0,20],[1;2;3]);利用命令plot(T,Y(:,1))可得积分曲线.例8.6求解Vander

Pol方程x''+(x2-1)x'+x=0

解：令y1=x,y2=y1'=x',则方程化为方程组

y1'

=y2,y2'

=(1-y12)y2-y1建立函数M文件vdp：function

dy=vdp(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=(1-y(1)^2)*y(2)-y(1);取初始条件y(0)=[0;1],Ts=[0,20],利用ODE45命令求解如下Ts=[0,20];y0=[0;1];[T,Y]=ode45('vdp',Ts,y0);积分曲线

由y1=x可知，原方程的解是输出矩阵Y的第一列Vander

Pol方程x''+(x2-1)x'+x=0满足