微分方程的基础知识与练习

微分方程的基本知识与练习（一）微分方程基本概念：一方面通过一种具体的问题来给出微分方程的基本概念。（1）一条曲线通过点（1，2），且在该曲线上任一点M（x,y）处的切线的斜率为2x，求这条曲线的方程。解 设曲线方程为.由导数的几何意义可知函数满足 （1）同步还满足如下条件：时， （2）把（1）式两端积分，得 即 （3）其中C是任意常数。把条件（2）代入（3）式，得， 由此解出C并代入（3）式，得到所求曲线方程： （4）（2）列车在水平直线路上以20的速度行驶；当制动时列车获得加速度.问开始制动后多少时间列车才干停住，以及列车在这段时间里行驶了多少路程？解 设列车开始制动后t秒时行驶了s米。根据题意，反映制动阶段列车运动规律的函数满足： （5）此外，还满足条件：时， （6）(5)式两端积分一次得： （7）再积分一次得 （8）其中都是任意常数。把条件“时”和“时”分别代入（）式和（）式，得把的值代入（7）及（8）式得 （9） （10）在（9）式中令，得到列车从开始制动到完全停止所需的时间：。再把代入（10）式，得到列车在制动阶段行驶的路程上述两个例子中的关系式（1）和（5），（6）都具有未知函数的导数，它们都是微分方程。1微分方程的概念一般地，凡具有未知函数、未知函数的导数及自变量的方程，叫做微分方程。未知函数是一元函数的方程叫做常微分方程；未知函数是多元函数的方程，叫做偏微分方程。我们只研究常微分方程。微分方程中所浮现的未知函数的最高阶导数的阶数，叫做微分方程的阶。例如，方程（1）是一阶微分方程；方程（5）是二阶微分方程方程。又如，方程是四阶微分方程。一般地，阶微分方程的形式是 （11）其中F是个变量的函数。这里必须指出，在方程（11）中，是必须浮现的，而等变量则可以不浮现。例如阶微分方程中，除外，其她变量都没有浮现。由前面的例子我们看到，在研究某些实际问题时，一方面要建立微分方程，然后找出满足微分方程的函数，就是说，找出这样的函数 ，把这函数代入微分方程能使该方程成为恒等式。这个函数就叫做该微分方程的解。例如，函数（3）和（4）都是微分方程（1）的解；函数（8）和（10）都是微分方程（5）的解。如果微分方程的解中具有任意常数，且任意常数的个数与微分方程的阶数相似，这样的解叫做微分方程的通解。例如，函数（3）是方程（1）的解，它具有一种任意常数，而方程（1）是一阶的，因此函数（3）是方程（1）的通解。又如，函数（8）是方程的解，它具有两个任意常数，而方程（5）是二阶的，因此函数（8）是方程（5）的通解。由于通解中具有任意常数，因此它还不能完全拟定地反映某一客观事物的规律性，必须拟定这些常数的值。为此，要根据问题的实际状况提出拟定这些常数的条件。例如，例1中的条件（2），例2中的条件（6），便是这样的条件。设微分方程中的未知函数为，如果微分方程是一阶的，一般用来拟定任意常数的条件是时，或写成 其中，都是给定的值；如果微分方程是二阶的，一般用来拟定任意常数的条件是：时，或写成 ，其中，和都是给定的值。上述条件叫做初始条件。拟定了通解中的任意常数后来，就得到了微分方程的特解。例如（4）式是方程（1）满足条件（2）的特解；（10）式是方程（5）满足条件（6）的特解。求微分方程满足初始条件的特解这样一种问题，叫做一阶微分方程的初值问题，记作 （13）二阶微分方程的初值问题是3、 例题例1 验证：函数 （14）是微分方程 （15）的解。解 求出所给函数（14）的导数 把及的体现式代入方程（15）得+函数（14）及其导数代入方程（15）后成为一种恒等式，因此函数（14）是微分方程（15）的解。用程序来实现： syms k t C1 C2; x=C1*cos(k*t)+C2*sin(k*t); diff(x,t,2)+k2*xans =k2*(C1*cos(k*t) + C2*sin(k*t) - C1*k2*cos(k*t) - C2*k2*sin(k*t) simple(ans)（二）微分方程的解一、几种会用到的函数：1、solve函数：Matlab中solve函数重要是用来求解线性方程组的解析解或者精确解。solve函数的语法定义重要有如下四种：solve(eq)solve(eq, var)solve(eq1,eq2, , eqn)g = solve(eq1, eq2, , eqn, var1, var2, , varn)eq代表字符串形式的方程，var代表的是变量。例1：解方程程序是：syms a b c x; solve(a*x2+b*x+c) （ 也可写成solve(a*x2+b*x+c=0) ）当没有指定变量的时候，matlab默认求解的是有关x的解，求解的成果为：ans = -(b + (b2 - 4*a*c)(1/2)/(2*a) -(b - (b2 - 4*a*c)(1/2)/(2*a)d当指定变量为b的时候：solve(a*x2+b*x+c,b)求解的成果为： ans =-(a*x2 + c)/xs = -(a*x2 + c)/x例2：对于方程组的状况S=solve(x+y=1,x-11*y=5);S.xS.y S=S.x,S.y(这里或者写成x=S.x y=S.y) 如果解得是一种方程组，并且采用了形如a,b=solve(a+b=1, 2a-b=4ab) 的格式，那么，在MATLAB Ra中没问题，可以保证输出的a，b就等于相应的解，但是在Rb等早先版本中不能保证输出的顺序就是你声明变量时的顺序。因此最佳采用g=solve(a+b=1, 2a-b=4ab)这种单输出格式，这样输出的是一种构造体，g.a和g.b就是相应的解。S = 4/3, -1/3一、 微分方程的解析解格式：dsolve(方程1, 方程2,方程n, 初始条件, 自变量) 记号: 在体现微分方程时，用字母D表达求微分，D2y、D3y 等表达求高阶微分.任何D后所跟的字母为因变量，自变量可以指定或由系统规则选定为确省，默认自变量是t例如，微分方程 应体现为：D2y=0.例1：求解微分方程，并加以验证求解本问题的Matlab 程序为：syms x y %line1y=dsolve(Dy+2*x*y=x*exp(-x2),x) %line2diff(y,x)+2*x*y-x*exp(-x2) %line3simplify(diff(y,x)+2*x*y-x*exp(-x2) %line4阐明：(1) 行line1是用命令定义x,y为符号变量这里可以不写，但为保证对的性，建议写上；(2) 行line2是用命令求出的微分方程的解：1/2*exp(-x2)*x2+exp(-x2)*C1(3) 行line3使用所求得的解这里是将解代入原微分方程，成果应当为0，但这里给出：-x3*exp(-x2)-2*x*exp(-x2)*C1+2*x*(1/2*exp(-x2)*x2+exp(-x2)*C1)(4) 行line4 用 simplify() （simple()）函数对上式进行化简，成果为 0， 表白的确是微分方程的解例2：先求微分方程的通解，再求在初始条件下的特解，并画出特解函数的图形求解本问题的 Matlab 程序为：syms x yy=dsolve(x*Dy+y-exp(x)=0, x)成果y =(exp(x)+C1)/x求特解两个措施1.y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1), x)成果y = (exp(x)+exp(1)/x2. C1= solve(2*exp(1)=exp(1)+C1,C1)成果C1 =exp(1)y =(exp(x)+exp(-x2)成果(exp(x)+exp(1)/xezplot(y) 例3：求微分方程组在初始条件下的特解，并画出解函数的图形求解本问题的 Matlab 程序为：syms x y ta=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=0,x(0)=1,y(0)=0,t); x=a.xy=a.ysimple(x);simple(y);ezplot(x,y,0,1.3);axis auto %坐标刻度选默认值例4 先求微分方程的通解，再求微分方程的特解. 程序是：dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x)ans = (3*sin(5*x)/exp(2*x)例5 求微分方程组的通解. 程序是：A=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z,t); x=A.xy=A.y z=A.z练习：求下列微分方程的解析解（1）程序：syms a b y;dsolve(Dy=a*y+b,x)ans =-(b - C2*exp(a*x)/a（2）程序： syms x y; dsolve(D2y=sin(2*x)-y,y(0)=0,Dy(0)=1,x)ans =(5*sin(x)/3 - sin(2*x)/3（3）程序：syms f g x;g=dsolve(Df=f+g,Dg=g-f,Df(0)=1,Dg(0)=1,t)g.g =exp(t)*sin(t)g.f =exp(t)*cos(x)（4）求微分方程的通解程序及答案：dsolve(x2-1)*Dy+2*x*y-sin(x)=0,x) ans = (-cos(x)+C1)/(x2-1)（5）求微分方程的通解程序及答案：dsolve(D2y-2*Dy+5*y=exp(x)*sin(x),x) ans = exp(x)*sin(2*x)*C2+exp(x)*cos(2*x)*C1+1/3*exp(x)*sin(x)（6）求微分方程组 在初始条件下的特解，并画出解函数的图形程序及答案：x,y= dsolve(Dx+x+y=0,Dy+x-y=0,x(0)=1,y(0)=0,t);或S= dsolve(Dx+x+y=0,Dy+x-y=0,x(0)=1,y(0)=1,t);x=S.xy=S.yezplot(x,y)解下列微分方程：二、 微分方程的数值解在生产和科研中所解决的微分方程往往很复杂且大多得不出一般解。而在事实上对初值问题，一般是规定得到解在若干个点上满足规定精确度的近似值，或者得到一种满足精确度规定的便于计算的体现式。因此，研究常微分方程的数值解法是十分必要的。一阶常微分方程与高阶微分方程可以互化，已给一种阶方程设，可将上式化为一阶方程组t，y=solver（f,ts,y0）ode45 ode23 ode113ode15sode23s由待解方程写成的m-文献名ts=t0，tf，t0、tf为自变量的初值和终值函数的初值自变量值函数值ode23：组合的2/3阶龙格-库塔-芬尔格算法ode45：运用组合的4/5阶龙格-库塔-芬尔格算法等等。（t,y=ode45(yprime,t0,tf,y0) ,yprime是用以表达f(t,y)的M文献名，t0表达自变量的初始值，tf表达自变量的终值，y0表达初始向量值。输出向量t表达节点(t0;t1;,tn),输出矩阵y表达数值解，每一列相应y的一种分量。ode45是最常用的求解微分方程数值解的命令，对于刚性方程组不适宜采用。ode23与ode45类似，只是精度低某些。ode12s用来求解刚性方程组，是用格式同ode45。求解器ODE类型特点阐明ode45非刚性一步算法，4,5阶Runge-Kutta措施累积截断误差大部分场合的首选算法ode23非刚性一步算法，2,3阶Runge-Kutta措施累积截断误差使用于精度较低的情形ode113非刚性多步法，Adams算法，高下精度均可达到计算时间比ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法，Gears反向数值积分，精度中档若ode45失效时，可尝试使用ode23s刚性一步法，2阶Rosebrock算法，低精度。当精度较低时，计算时间比ode15s短注意:1、在解n个未知函数的方程组时，y0和y均为n维向量，m-文献中的待解方程组应以y的分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.2,例1：求解微分方程初值问题的数值解，求解范畴为区间0, 15程序是：先建立M文献(函数M 文献的格式：function 返回变量=函数名（自变量，因变量）注释阐明语句段程序语句段）可以编这样的程序：function dy=a1(x,y)dy=-2*y+2*x2+2*x;x,y=ode45(a1,0,15,1)plot(x,y,ro)（也可以编下面的程序function dy=ac(x,y) dy=zeros(1,1); dy=-2*y+2*x2+2*x;再编写程序x,y=ode45(ac,0,15,1)plot(x,y,o)）例令 y1=x，y2=y1化为一阶微分方程组：程序1：先建立M文献function dy=a2(t,y)dy=y(2);1000*(1-y(1)2)*y(2)-y(1)或function dy=a3(t,y)dy(1)=y(2);dy(2)=1000*(1-(y(1)2)*y(2)-y(1);dy=dy(1);dy(2);输入命令t,y=ode15s(a2,0,1000,2,0)plot(t,y(:,1)程序2：先建立M文献function dy=xx(t,y)dy=zeros(2,1)dy(1)=y(2);dy(2)=1000*(1-y(1)2)*y(2)-y(1);再编写程序t,y=ode15s(xx,0,1000,2,0); plot(t,y(:,1),r*)选定a=0.3，b=2,c=3 初值x(0)=0,y(0)=0,z(0)=0化为一阶微分方程组：令程序1:先建立M文献function dy=ros(t,y)dy(1)=-y(2)-y(3)dy(2)=y(1)+0.3*y(2)dy(3)=2+(y(1)-3)*y(3)dy=dy(1);dy(2);dy(3)再编写程序：t,y=ode45(ros,0,10,0,0,0); plot(t,y(:,1),r,t,y(:,2),g,t,y(:,3),y)程序2：先建立M文献function dy=ros(t,y)a=0.3;b=2;c=3;dy=zeros(3,1);dy(1)=-y(2)-y(3);dy(2)=y(1)+a*y(2);dy(3)=b+(y(1)-c)*y(3);再编写程序t,y=ode45(ros,0,10,0,0,0); plot(t,y(:,1),r,t,y(:,2),g,t,y(:,3),y)例4、求解微分方程先求解析解，再求数值解，并画图进行比较。先求数值解程序1：先建立M文献function dy=jie(t,y)dy=-y+t+1;再编程t,y=ode45(jie,0,10,1); plot(t,y,r*)hold on t=0:0.5:10; dsolve(Dy=-y+t+1,y(0)=1)程序2：先建立M文献function dy=jie(t,y)dy=zeros(1,1);dy=-y+t+1;再编程t,y=ode45(jie,0,10,1); plot(t,y,r*)hold on t=0:0.5:10; dsolve(Dy=-y+t+1,y(0)=1)y=t+1./exp(t); plot(t,y) xlabel(t),ylabel(y)例6、分别用 ode23、ode45 求 微分方程初值问题的数值解(近似解)，求解区间为运用画图来比较两种求解器之间的差别解：令，则方程变为：解：先建立M文献function dy=a(t,y)dy(1)=-y(1)-y(2);dy(2)=-y(1)+y(2);dy=dy(1);dy(2);再编程序t,y=ode23(a,0,2,1,0) plot(t,y(:,1),r*,t,y(:,2),g+) hold on; t,y=ode45(a,0,2,1,0); plot(t,y(:,1),b,t,y(:,2),y)例 求解微分方程先求解析解，再求数值解，并进行比较。由clear; s=dsolve(Dy=-y+t+1,y(0)=1,t)simplify(s)可得解析解为。下面再求其数值解，先编写M文献fun8.m%M函数fun8.mfunction f=fun8(t,y)f=-y+t+1;(红色部分用来求数值解，与解析解无关）再用命令clear; close; t=0:0.1:1;y=t+exp(-t); plot(t,y); %化解析解的图形hold on; %保存已经画好的图形，如果下面再画图，两个图形和并在一起t,y=ode45(fun8,0,1,1); plot(t,y,ro); %画数值解图形，用红色小圈画xlabel(t),ylabel(y)成果见图7.1图16.6.1 解析解与数值解由图16.6.1可见，解析解和数值解吻合得较好。例3 求方程的数值解.不妨取.则上面方程可化为先看看有无解析解.运营MATLAB代码s=dsolve(D2y=9.8*sin(y),y(0)=15,Dy(0)=0,t) %该命令输入后显示无解析解知原方程没有解析解.下面求数值解.令可将原方程化为如下方程组建立M文献a5.m如下% function dy=a5(t,y)dy(1)=y(2);dy(2)=9.8*sin(y(1);dy=dy(1);dy(2)运营MATLAB代码clear; close; t,y=ode45(a5,0,10,15,0); plot(t,y(:,1); %画随时间变化图,y(:2)则表达的值xlabel(t),ylabel(y1)成果见图7.2图7.2 数值解练习：1求方程的解析解和数值解,并进行比较程序：先建立M 文献function dy=a2(x,y)dy=y(2);(2*x*y(2)/(1+x2)再输入程序x,y=ode45(a2,0,10,1,3)plot(x,y(:,1),r+)hold ony=dsolve(1+x2)*D2y=2*x*Dy,y(0)=1,Dy(0)=3,x)ezplot(y)或者hold onx=0:0.1:10;y=1+3.*x+x.3;plot(x,y,g*)2分别用ode45和ode15s求解Van-del-Pol方程的数值解和解析解,并进行比较.（该题没有解析解）设x=x(1)，程序：先建立M文献function dx=a3(t,x)dx=x(2);1000*(1-x(1)2)*x(2)+x(1)再输入程序t,x=ode45(a3,0,1,0,1);plot(t,x(:,1),r*) hold ont,x=ode15s(a3,0,1,0,1);plot(t,x(:,1),bo)x=dsolve(D2x-1000*(1-x2)*Dx-x=0,x(0)=0,Dx(0)=1,t) %此命令输出成果显示没有解析解