科学计算方法20(常微分方程数值解)

1、引例1. 人口模型,引例1. 人口模型,引例1. 人口模型,引例1. 人口模型,差分方法,(离散化),难,易,(连续化),插值,数值方法,前向差分公式,后向差分公式,中心差分公式,例1. 用Euler法求初值问题的数值解。,解: 步长h=0.1, xn= nh (n = 0, 1, 10), 将xn 代入方程,Euler公式: yn+1 = yn + h( yn 2xn /yn) (n = 0, 1, ,10),并用前向差分格式代替其中的导数项,方程中含有导数项y, 这是微分方程的本质特征, 也是它难以求解的症结所在。常见解决思路通常为数值微分和数值积分。,将xn 代入方程,并用前向差分格式代

2、替其中的导数项,若用y(xn)的近似值yn代入上式, 并记所得结果为yn+1,Euler公式:,例1. 用Euler法求初值问题的数值解。,解: 步长h=0.1, xn= nh (n = 0, 1, 10),Euler公式: yn+1 = yn + 0.1( yn 2xn /yn) (n = 0, 1, ,10),clear; f = inline(y-2*x/y,x,y); a = 0; b = 1;n=100; h =(b-a)/n; x=a:h:b; y(1) = 1; for i = 1 : n y(i+1) = y(i) + h*f(x(i),y(i); end plot(x,y,m

3、s-); hold on, y=(1+2*x).(1/2); plot(x,y),方程中含有导数项y, 这是微分方程的本质特征, 也是它难以求解的症结所在。,将xn+1 代入方程,并用后向差分格式代替其中的导数项,若用y(xn)的近似值yn代入上式, 并记所得结果为yn+1,隐式Euler公式:,y = f (x, y),左矩形积分公式,Euler公式：,方程中含有导数项y, 这是微分方程的本质特征, 也是微分方程难以求解的症结所在。常见解决思路通常为数值微分和数值积分。,若用y(xn)的近似值yn代入上式, 并记所得结果为yn+1,y = f (x, y),右矩形积分公式,若用y(xn)的近

4、似值yn代入上式, 并记所得结果为yn+1,隐式Euler公式:,y = f (x, y),若用y(xn)的近似值yn代入上式, 并记所得结果为yn+1,梯形积分公式:,梯形公式:,显式格式与隐式格式各有利弊: 显式格式的计算量小, 但稳定性较差; 与此相反, 隐式格式的稳定性好, 但需要迭代求解, 计算量比较大。,预报,预报-校正方法(改进Euler公式):,综合使用这两种方法, 先用Euler公式求得一个初步的近似值, 称为预报值; 预报值的精度不高, 代入右端的yn+1计算校正值。,校正,梯形公式：,改写如下 K1 = f(xn , yn) , K2 = f( xn+1 , yn+ hK

5、1) ,改进Euler公式:,例1. 用改进Euler法求初值问题的数值解。,解: 步长h=0.1, xn= nh (n = 0, 1, 10),clear; f = inline(y-2*x/y,x,y); a = 0; b = 1;n=10; h =(b-a)/n; x=a:h:b; y(1) = 1; for i = 1 : n K1=f(x(i),y(i); K2=f(x(i+1),y(i)+h*K1); y(i+1) = y(i) + (h/2)*(K1+K2); end hold on, plot(x,y,gs-); hold on, y=(1+2*x).(1/2); plot(x

6、,y),计算方法：算法设计及其Matlab实现 王能超编,构造微分方程数值格式就是估计斜率K=f(,y()。,根据微分中值定理,Euler公式: yn+1 = yn+ hf (xn, yn),隐式Euler公式: yn+1 = yn+ hf (xn+1, yn+1),改写如下 K1 = f(xn , yn) , K2 = f( xn+1 , yn+ hK1) ,改进Euler公式:,对区间xn, xn+1内xn和xn+p=xn+ph (0p 1)斜率值K1和K2加权平均,问题是如何生成y(xn+p)?,设yn已知,则Euler公式预报 y(xn+p),设计格式如下: K1 = f(xn , y

7、n) , K2 = f( xn+ph, yn+ phK1) ,我记得我的朋友冯诺依曼曾经说过,用四个参数我可以拟合出一头大象,而用五个参数我可以让它的鼻子摆动。 费米,其中,例1. 用四级Range-Kutta公式求初值问题的数值解。,解:,K1=f(xn, yn), K2=f(xn+0.5h, yn+0.5hK1) K3=f(xn+0.5h, yn+0.5hK2), K4=f(xn+h, yn+hK3),K1=yn 2xn/yn, K2= (yn+0.5hK1) 2(xn+0.5h)/(yn+0.5hK1) K3= (yn+0.5hK2) 2(xn+0.5h)/(yn+0.5hK2) K4=

8、 (yn+hK3) 2(xn+h)/(yn+hK3),例1. 用四级Range-Kutta公式求初值问题的数值解。,解： 步长h=0.2, xn= nh (n = 0, 1, 5),clear f = inline(y-2*x/y,x,y); a = 0; b = 1; n = 5; h = (b-a)/n; x = a : h : b; y(1) = 1; for i = 1 : n K1 = f(x(i), y(i); K2 = f(x(i)+h/2, y(i)+K1*h/2); K3 = f(x(i)+h/2, y(i)+K2*h/2); K4 = f(x(i)+h, y(i)+K3*h

9、); y(i+1) = y(i) + h*(K1+2*K2+2*K3+K4)/6; end plot(x,y,rs-); %hold on, %y=(1+2*x).(1/2); %plot(x,y),Runge-Kutta 方法,Matlab 解初值问题函数,用 Maltab自带函数 解初值问题,求解析解：dsolve,求数值解： ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb,dsolve 求解析解,dsolve 的使用,y=dsolve(eq1,eq2, . ,cond1,cond2, . ,v),其中 y 为输出的解, eq1、eq2、

10、. 为微分方程, cond1、cond2、. 为初值条件, v 为自变量。,例 ：求微分方程 的通解, 并验证。, y=dsolve(Dy+2*x*y=x*exp(-x2),x), syms x; diff(y)+2*x*y - x*exp(-x2),dsolve 的使用,几点说明,如果省略初值条件, 则表示求通解;,如果省略自变量, 则默认自变量为 t,dsolve(Dy=2*x,x); % dy/dx = 2x dsolve(Dy=2*x); % dy/dt = 2x,微分方程中用 D 表示对 自变量 的导数, 如：,Dy y； D2y y； D3y y,数值求解,T,Y = solver

11、(odefun,tspan,y0),其中 y0 为初值条件, tspan为求解区间; Matlab在数值求解时自动对求解区间进行分割, T (列向量) 中返回的是分割点的值(自变量), Y (数组) 中返回的是这些分割点上的近似解, 其列数等于因变量的个数。 solver 为Matlab的ODE求解器(可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb),没有一种算法可以有效地解决所有的 ODE 问题, 因此MATLAB 提供了多种ODE求解器, 对于不同的ODE, 可以调用不同的求解器。,Matlab的ODE求解器,Reference: T

12、he Matlab ODE suite, SIAM J. Scientific Computing,假如一个被求取的解变化缓慢，但隔不久又会突然快速变化，于是数值解算法必须取很小的步长才能获得满意的结果，那么这种问题就是刚性问题。刚性问题就是效率问题，假如我们不关心计算耗费时间的长短，就不必理会刚性。非刚性算法可以求解刚性问题，只不过需要花费很长的时间。,参数说明,odefun 为显式常微分方程, 可以用命令 inline 定义, 或在函数文件中定义, 然后通过函数句柄调用。,fun=inline(-2*y+2*x2+2*x,x,y); x,y=ode23(fun,0,0.5,1);,T,Y = solver(odefun,tspan,y0),引例2. Lorenz模型 (lorenzgui),根据大气运动的规律,Lorenz建立了简化的数学模型。Lorenz经过研究发现,天气预测具有对初始条件的敏感依赖性, 即初始条件最微小的差异都会导致天气的行为无法准确预测。Lorenz结论說天气的长期预报是不可能的。1979年12月, Lorenz在华盛顿