求导积分与微分方程数值解(第2次课.ppt

第三讲 微分方程数值解,内容：本讲首先以单摆微分方程求解过程引入， 介绍基于MATLAB的微分方程求解函数， 然后重点讲解微分方程(组) 图形图像解法， 最后简介 Differential Equation Editor 目的：掌握微分方程数值解的一般思路和方法 要求：能够处理应用类型微分方程数值解问题 掌握Malthus人口模型参数计算 (本讲实验题目) 掌握常用求解函数 dsolve ode23 ode45 ode15s 掌握图形图像求解方法 斜率场 / 相平面 / 等值线 了解基于Simulink的微分方程数值仿真工具 DEE,大多数微分方程无法求解析解?,微分方程是研究函数变化规律的有力工具，在科技、工程、经济管理、及生态、环境、人口、交通等各个领域有着广泛的应用。 建立微分方程可以依据物理的、或其他原理和规律建立的平衡关系，但是！更重要的问题是如何求解这些微分方程(组) 部分微分方程可以求得解析解，但是绝大多数的非线性、变系数微分方程或“难以求解”或“求不出解”，所以对于实际问题，研究微分方程的数值解具有重要意义！,ODE vs PDE .,我们所熟悉的微分方程？,微分方程初值问题的最简单形式：,简单定义：含有导数的方程就称为微分方程 一般形式： 解析解：求得具体解析式y=f(x) (早先学过的.) 数值解：求得系列散点xi对应的近似值yi(表格法) 图像解：用图像表示解曲线(有何优势？),表示函数的三种方法？,单摆微分方程求解：建立方程,引例：单摆微分方程求解，由已知运动规律建立微分方程，在求解过程中采用了两种方法，方法一近似简化，方法二求数值解 (实验室练习),由牛顿第二定律建立方程：,我们要找到符合条件的theta与t的函数关系，但是此2阶非线性微分方程并不易求得解析解除非,单摆微分方程求解：求近似解,除非一：简化方程求其近似解： 取x0=10即x0= 0.1745，在此弧度范围内sin 所以原微分方程可以简化为：,此线性常系数微分方程可以直接用dsolve函数求得：,dsolve(D2theta+g/l*theta=0,theta(0)=a0,Dtheta(0)=0,t) 其解析解为： a0*cos(g/l)(1/2)*t) 这个解可以作为原方程的近似解,简化,单摆微分方程求解：求数值解,除非二：求其数值解，也就是部分点xi对应的近似值yi 我们采用ode23函数求解，所以首先改写方程：,改写,由简化方程求得近似解为：,代入g=9.8 l=25，得到周期 T 10s，下面考察ts=0到tf=10内若干点处的近似值i，即所谓的数值解,单摆微分方程求解：求数值解,首先建立被调函数danbai.m function xdot=danbai(t,x) g=9.8;l=25; xdot(1)=x(2); xdot(2)=-g/l*sin(x(1);xdot=xdot;,然后是主调指令，也可写成主调文件loaddanbai.m warning off ts=0; tf=10; a0=0.1745; cond0=a0,0; %初始化变量 t,x=ode23(danbai,ts,tf,cond0); %调用ode23函数求解 g=9.8; l=25; w=sqrt(g/l); y=a0*cos(w*t); %近似解 t,x(:,1),y %输出t对应的数值解和近似解 subplot(1,2,2); stem(t,x(:,1),ro); title(数值解) subplot(1,2,1); hold on; stem(t,y, bp); plot(t,y, b-); title(近似解),用dsolve函数求解微分方程,MATLAB求解微分方程解析解的函数dsolve Symbolic solution of ordinary differential equations . Syntax r = dsolve(eq1,eq2,.,cond1,cond2,.,v) 题例1：p49-4.4.1/ex1，ex2 dsolve(Dy=1+y2) dsolve(Dtheta=1+theta2,theta(0)=1,xi) dsolve(x2*D2y+x*Dy+(x2-(1/2)2)*y=0,y(pi/2)=2,Dy(pi/2)=-2/pi,x) pretty(ans) 提示：一些需要注意的细节,用dsolve函数求解微分方程,MATLAB求解微分方程组解析解的函数dsolve 题例2：求解 dsolve(Dy+2*x*y=x*exp(-x2),x) 题例3：求解 dsolve(x2-1)*Dy+2*x*y-cos(x)=0,y(0)=1,x) 题例4：求解 dsolve(D2y+3*Dy+exp(x)=0,x),用dsolve函数求解微分方程组,MATLAB求解微分方程组解析解的函数dsolve 题例5： p50-4.4.1/ex3，ex4 f,g=dsolve(Df=3*f+4*g,Dg=-4*f+3*g) f,g=dsolve(Df=3*f+4*g,Dg=-4*f+3*g,f(0)=0,g(0)=1,x) 下面的指令有否区别？ dsolve(Dy=x*sin(x)/cos(y) dsolve(Dy=x*sin(x)/cos(y),x) 提示： 用dsolve求解存在解析解的微分方程相当方便，在“只要结果，不求过程”的场合，节约了大量时间。,课本引例: Malthus人口模型计算,Malthus认为单位时间内人口净增长率为常数：,d=1790:10:1900; x=3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76; t=(d-1790)/10; y=log(x); p=polyfit(t,y,1); y=poly2str(p,t) r=p(1), a=p(2); x0=exp(a) plot(t,x,r+); hold on; t=min(t):0.01:max(t); x=x0*exp(r*t); plot(t,x,b-);,微分方程数值的相关方法,求微分方程数值解的方法很多，比如：欧拉法，龙格-库塔法等。Runge-Kutta法本质上是间接使用Taylor级数的一种技术，有兴趣的同学可以参考微分方程数值解相关专著。 首选ode45 ，次选ode15s，低精度可选ode23,基于R-K算法的数值解函数ode23,MATLAB求解微分方程数值解的函数ode23 Solve initial value problems for ordinary differential equations (ODEs).Syntax(ode23): T,Y = ode23(odefun,tspan,y0) T,Y = ode23(odefun,tspan,y0,options) 题例6：p51-ex1,基于R-K算法的数值解函数ode45,MATLAB求解微分方程数值解的函数ode45 T,Y = ode45(odefun,tspan,y0) T,Y = ode45(odefun,tspan,y0,options) 题例7：p56-ex,基于R-K算法的数值解函数ode45,MATLAB求解微分方程数值解的函数ode45 T,Y = ode45(odefun,tspan,y0) T,Y = ode45(odefun,tspan,y0,options) 题例8：Lorenz微分方程模型 loadlrz(0,50,0,0,0.1,8/3,10,28) loadlrz(0,30,0.1,0.1,0.21,3.1,23,8) loadlrz(0,100,0.1,0.1,0.21,3.92,11,28),解析法/数值法的局限.,解析法能得到精确解，能绘制图形，但适应面窄; 数值法能得到离散点的近似值，但无法兼顾整体; 图像法在数值法的基础上更好地展示全局和趋势. 方法概述 ： 斜率场法是依据y=f(x,y)得到平面上一些点的斜率值，然后过这些点引出自该点出发的短直线，通过观察趋势，了解解曲线的分布和性态。 题例9：p45-ex1,图解法-斜率场1,用斜率场法求解微分方程：y=sinx siny syms x y; fun=sin(x)*sin(y); hx=16/40; hy=16/40; x0=-8; y0=-8; hold on for i=1:40 x=x0+(i-1)*hx; for j=1:40 y=y0+(j-1)*hy; k=eval(fun); if abs(hx*k)hy plot(x,x+hy/k*2/3,y,y+hy*2/3) else plot(x,x+hx*2/3,y,y+hx*k*2/3) end end end title(dy/dx=sinx*siny); xlabel(x); ylabel(y);,图解法-斜率场2,图解法-相平面轨迹1, 方法概述 ： 相平面轨迹法是依据不同的初值条件，先用数值解法求出各自对应的数值解（x(t),y(t)），最后用plot描点绘图，通过观察趋势，了解解曲线的分布和性态。 题例10：p53-ex1 先用数值解法求出若干初值条件下的(x(t),y(t) %tbp53.m function dequ=tbp53(t,x) dequ=2*x(1)-1*x(1)*x(2);-1*x(2)+1*x(1)*x(2);,图解法-相平面轨迹2,%loadtbp53.m hold on for i=1:7 tspan=0:0.01:5; cond0=1,0.1+0.2*i; t,x=ode45(tbp53,tspan,cond0); plot(x(:,1),x(:,2) end axis(0 4 0 4); xlabel(x); ylabel(y);,图解法-等值线, 方法概述 : 等值线隶属于相平面轨迹法，先求出通解(dsolve？)，再针对通解中的常数，每一个常数定值都对应着一条等值线，用contour函数根据三维数据绘出等值线即可，通过观察趋势，了解解曲线的分布和性态。 题例11：p55- x,y=meshgrid(0:.1:4, 0:.1:4); z=2*log(y)-y+log(x)-x; contour(x,y,z,20); xlabel(x);ylabel(y),简介微分方程编辑器 DEE 1,基于Simulink的微分方程数值仿真工具 DEE 题例12： x,y=dsolve(Dx=x+y*exp(-t),Dy= -x*y+cos(t),x(0)=0,y(0)=1,t) %Explicit solution could not be found DEE 1 x(1)+x(2)*exp(-u(1) -x(1)*x(2)+cos(u(1) 0 1 x(1) x(2),简介微分方程编辑器 D