《求微分方程的解》PPT课件.ppt

1,实验四 求微分方程的解,数学实验,2,自牛顿发明微积分以来，微分方程在描述事物运动规律上已发挥了重要的作用。实际应用问题通过数学建模所得到的方程，绝大多数是微分方程。,由于实际应用的需要，人们必须求解微分方程。然而能够求得解析解的微分方程十分有限，绝大多数微分方程需要利用数值方法来近似求解。,本实验主要研究如何用 Matlab 来计算微分方程（组）的数值解，并重点介绍一个求解微分方程的基本数值解法Euler折线法。,问题背景和实验目的,3,考虑一维经典初值问题,基本思想：用差商代替微商,根据 Talyor 公式，y(x) 在点 xk 处有,Euler 折线法,4,初值问题的Euler折线法,具体步骤：,等距剖分：,步长：,分割求解区间,分割求解区间，差商代替微商，解代数方程,为分割点,5,Euler 折线法举例,例：用 Euler 法解初值问题,取步长 h = (2 - 0)/n = 2/n，得差分方程,当 h=0.4，即 n=5 时，Matlab 源程序见 fuluA.m,解：,6,Euler 折线法源程序,clear f=sym(y+2*x/y2); a=0; b=2; h=0.4; n=(b-a)/h+1; % n=(b-a)/h; x=0; y=1; szj=x,y; for i=1:n-1 % i=1:n y=y+h*subs(f,x,y,x,y); x=x+h; szj=szj;x,y; end szj plot(szj(:,1),szj(:,2),or-),7,Euler折线法举例（续）,解析解：,解析解,近似解,y=1/3*(-18-54*x+45*exp(3*x)(1/3),8,Runge-Kutta 方法,为了减小误差，可采用以下方法：,让步长 h 取得更小一些；,改用具有较高精度的数值方法：,龙格-库塔方法,Runge-Kutta (龙格-库塔) 方法,是一类求解常微分方程的数值方法,有多种不同的迭代格式,9,Runge-Kutta 方法,用得较多的是 四阶R-K方法,其中,10,四阶 R-K 方法源程序,clear; f=sym(y+2*x/y2); a=0; b=2; h=0.4; n=(b-a)/h+1; % n=(b-a)/h; x=0; y=1; szj=x,y; for i=1:n-1 % i=1:n L1=subs(f,x,y,x,y); L2=subs(f,x,y,x+h/2,y+L1*h/2); L3=subs(f,x,y,x+h/2,y+L2*h/2); L4=subs(f,x,y,x+h,y+L3*h); y=y+h*(L1+2*L2+2*L3+L4)/6; x=x+h; szj=szj;x,y; end plot(szj(:,1),szj(:,2), dg-),11,Runge-Kutta 方法,12,Euler 法与 R-K法误差比较,13,Matlab 解初值问题,用 Maltab自带函数 解初值问题,求解析解：dsolve,求数值解： ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb,14,dsolve 求解析解,dsolve 的使用,y=dsolve(eq1,eq2, . ,cond1,cond2, . ,v),其中 y 为输出， eq1、eq2、.为微分方程，cond1、cond2、.为初值条件，v 为自变量。,例 1：求微分方程 的通解，并验证。, y=dsolve(Dy+2*x*y=x*exp(-x2),x), syms x; diff(y)+2*x*y - x*exp(-x2),15,dsolve 的使用,几点说明,如果省略初值条件，则表示求通解；,如果省略自变量，则默认自变量为 t,dsolve(Dy=2*x,x); dy/dx = 2x dsolve(Dy=2*x); dy/dt = 2x,若找不到解析解，则返回其积分形式。,微分方程中用 D 表示对 自变量 的导数，如：,Dy y； D2y y； D3y y,16,dsolve 举例,例 2：求微分方程 在初值条件 下的特解，并画出解函数的图形。, y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x) ezplot(y);,17,dsolve 举例,例3：求微分方程组 在初值条件 下的特解，并画出解函数的图形。,x,y=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=0, . x(0)=1, y(0)=0, t) ezplot(x,y,0,1.3);,注：解微分方程组时，如果所给的输出个数与方程个数相同，则方程组的解按词典顺序输出；如果只给一个输出，则输出的是一个包含解的结构（structure）类型的数据。,18,dsolve 举例,例：,x,y=dsolve(Dx+5*x=0,Dy-3*y=0, . x(0)=1, y(0)=1,t),r = dsolve(Dx+5*x=0,Dy-3*y=0, . x(0)=1, y(0)=1,t),这里返回的 r 是一个 结构类型 的数据,r.x %查看解函数 x(t) r.y %查看解函数 y(t),只有很少一部分微分方程（组）能求出解析解。 大部分微分方程（组）只能利用数值方法求数值解。,dsolve的输出个数只能为一个 或 与方程个数相等,19,数值求解,T,Y = solver(odefun,tspan,y0),其中 y0 为初值条件，tspan为求解区间；Matlab在数值求解时自动对求解区间进行分割，T (列向量) 中返回的是分割点的值(自变量)，Y (数组) 中返回的是这些分割点上的近似解，其列数等于因变量的个数。 solver 为Matlab的ODE求解器（可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）,没有一种算法可以有效地解决所有的 ODE 问题，因此MATLAB 提供了多种ODE求解器，对于不同的ODE，可以调用不同的求解器。,20,Matlab的ODE求解器,21,参数说明,odefun 为显式常微分方程，可以用命令 inline 定义，或在函数文件中定义，然后通过函数句柄调用。,fun=inline(-2*y+2*x2+2*x,x,y); x,y=ode23(fun,0,0.5,1);,注：也可以在 tspan 中指定对求解区间的分割，如：,x,y=ode23(fun,0:0.1:0.5,1); %此时 x=0:0.1:0.5,T,Y = solver(odefun,tspan,y0),22,数值求解举例,如果需求解的问题是高阶常微分方程，则需将其化为一阶常微分方程组，此时必须用函数文件来定义该常微分方程组。,令,23,数值求解举例,先编写函数文件 verderpol.m,function xprime=verderpol(t,x) global mu; xprime=x(2); mu*(1-x(1)2)*x(2) - x(1);,再编写脚本文件 vdpl.m，在命令窗口直接运行该文件,clear; global mu; mu=7; y0=1; 0; %t,x=ode45(verderpol, 0,40, y0); t,x=ode45(verderpol, 0,40, y0); plot(t,x(:,1);,24,求解微分方程小结,Matlab 函数,求解析解（通解或特解），用 dsolve 求数值解（特解），用 ode45、ode23 .,Matlab 编程,Euler 折线法 Runga-Kutta 方法,25,如何定义函数文件,T,Y = ode45(odefun,tspan,y0) T,Y = ode23(odefun,tspan,y0),当 是函数向量时呢？,那么odefun就是,26,如何定义函数文件,T,Y = ode45(odefun,tspan,y0) T,Y = ode23(odefun,tspan,y0),ode45、ode23 等函数可用于求解显式常微分方程 当 是向量函数时，所对应的方程即为微分方程组,odefun,27,举例说明,fun=inline(-2*y+2*x2+2*x,x,y); x,y=ode23(fun,0,0.5,1);,例：求初值问题 的数值解。,解法一：使用 inline 定义微分方程 odefun,odefun 为方程右端项 f(t,y) 可以用 inline 定义（只适合于单个方程的情形） 通过函数文件定义，然后用函数句柄调用（适合所有情形）,注：自变量必须在前面，因变量在后面！,28,举例说明（单个方程）,function dy = myfun1(x,y) dy = -2*y+2*x2+2*x;,解法二：通过函数文件定义微分方程 odefun,1、先编写函数文件 myfun1.m,clear; x,y=ode23(myfun1,0,0.5,1);,2、编写主文件 main1.m,或直接在 Matlab 命令窗口输入上面的语句。,29,举例说明（方程组）,解：此时只能通过函数文件定义微分方程 odefun,例：求 , , 的数值解。,function dy = myfun2(t,y) dy = zeros(3,1); % dy must be a column vector! dy(1) = y(2) * y(3); dy(2) = -y(1) * y(3); dy(3) = -0.51 * y(1) * y(2);,1、先编写函数文件 myfun2.m,clear; T,Y=ode45(myfun2,0,12,0,1,1);,2、编写主文件 main2.m,dy = y(2)*y(3); -y(1)*y(3); . -0.51*y(1)*y(2);,30,思考,function dy = myfun2(t,y) dy = zeros(3,1); dy(1) = y(2) * y(3); dy(2) = -y(1) * y(3); dy(3) = -0.51 * y(1) * y(2);,1、函数文件 myfun2.m,能不能写成下面形式？,function dy = myfun2(t,x,y,z) dy = zeros(3,1); dy(1) = y * z; dy(2) = -x * z; dy(3) = -0.51 * x * y;,X,31,说明,odefun,变量属性必须一一对应！,function dy = myfun2(t,y),如果是常微分方程组，y 就是列向量！,dy 必须是列向量，长度为方程组的个数，通常与y的长度相同！,函数中的输入参数和输出参数是形参，名字可以任意取，但必须满足上述条件。即输入参数有两个，第一个表示自变量，第二个是由因变量组成的列向量，输出参数必须是列向量。,32,例,function dy = myfun3(t,y) dy = zeros(2,1); dy(1) = y(2) + t; dy(2) = t - 2;,例：解初值问题： ， ，,function out = myfun3(t,y) out = y(2) + t; t 2;,function yprime = myfun3(x,y) yprime = y(2) + x; x 2;,clear; T,Y=ode45(myfun3,0,10,1,1);,2、主文件 main3.m,1、函数文件 myfun3.m,33,高阶常微分方程,高阶常微分方程,例： Van der Pol 初值问题,令 ，则原方程可化为,一阶常微分方程组,变量替换 化为,参数怎么处理？,用