微分方程数值解法(中科大MATLAB课程).pptVIP

微分方程数值解法 中国科技大学信息学院 陆伟 luwei@ustc.edu.cn 微分方程数值解法 Matlab中包含一套完整的微分方程数值解法的算法程序 ODE：常微分方程（初值） BVP：边界值问题的常微分方程 DDE：延迟微分方程 IDE： 隐微分方程 PDE：偏微分方程 其符号工具箱(Symbolic Math Toolbox)提供了求微分方程解析解的算法程序 微分方程数值解法 本节讲解Matlab中初值问题常微分方程的数值解法。 ODE: ordinary differential equations Initial Value Problems for ODEs 微分方程数值解法 微分方程数值解法介绍--- Euler 折线法 微分方程数值解法 基本思想：用差商代替微商 微分方程数值解法 微分方程数值解法 微分方程数值解法 a = 0; %[a b]为积分区间 b = 2; n = 5; %中间点数 h = (b-a)/(n-1); %步长 x = linspace(a,b,n); y = zeros(size(x)); y(1) = 1; for k = 1:n-1 y(k+1) = y(k) +h *(y(k)+x(k)/(y(k).^2)); end plot(x,y,o-); 微分方程数值解法 微分方程数值解法 步长减小，精度会提高。 微分方程数值解法 微分方程数值解法 微分方程数值解法 微分方程数值解法 微分方程数值解法 微分方程数值解法 微分方程数值解法 微分方程数值解法 微分方程数值解法 * * 考虑一维经典初值问题 方程 初值 求值区间 根据 Talyor 公式，y(x) 在点 xk 处有 具体步骤： 分割求解区间，差商代替微商，解代数方程 分割求解区间 等距剖分： 步长： 差商代替微商 得方程组： 例：用 Euler 法解初值问题 解： 取步长 h = (2 - 0)/n = 2/n，得差分方程 练习：编写matlab程序，用Euler法解上述微分方程 解析解： 为了减小误差，可采用以下方法： 让步长 h 取得更小一些； 改用具有较高精度的数值方法： Runge-Kutta (龙格-库塔) 方法 四阶R-K方法 其中 a = 0; b = 2; n = 50; h = (b-a)/(n-1); x = linspace(a,b,n); y = zeros(size(x)); y(1) = 1; f = @(u,v) u+v/u^2; for k = 1:n-1 L1 = f(y(k),x(k)); L2 = f(y(k)+h*L1/2, x(k)+h/2); L3 = f(y(k)+h*L2/2, x(k)+h/2); L4 = f(y(k)+h*L3, x(k)+h ); y(k+1) = y(k) +h *(L1+2*L2+2*L3+L4)/6; end 用 Maltab自带函数 解初值问题 求解析解：dsolve 求数值解：  ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb [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，可以调用不同的求解器。 使用于精度较低的情形 单步法；2，3 阶 R-K 方法；累计截断误差为 (△x)3 非刚性 ode23 计算时间比 ode45 短 多步法；Adams算法；高低精度均可到 10-3～10-6 非刚性 ode113 适度刚性情形 采用梯形算法 适度刚性 ode23t 若 ode45 失效时，可尝试使用 多步法；Gear’s 反向数值微分；精度中等 刚性 ode15s 当精度较低时，计算时间比 ode15s 短 单步法；2 阶Rosebrock 算法；低精度 刚性 ode23s 当精度较低时，计算时间比ode15s短 梯形算法；低精度 刚性 ode23tb 说明 特点 ODE类型 求解器 大部分场合的首选方法 单步法；4，5 阶 R-K 方法；累计截断误差为 (△x)3 非刚性 ode45 fun=inl