用matlab解微分方程初值问题,常微分方程初值问题的数值求解及MATLAB实现

科技信息SCIENCE & TECHNOLOGY INFORMATION2012 年 第 7 期 0 引言 在自然科学和经济生活的许多领域中,常常会遇到一阶常微分方程的初值问题如下：y'(t)=f(t，y)y(t0 )=y0 0 ， t0≤t (1) 这里 f 是凸区域 D 上的连续且满足 Lipschitz 条件的函数,可以证明常微分方程初值问题(1)具有唯一解. 1 常微分方程初值问题的数值解法 常微分方程初值问题的数值解法一般分为两大类，即单步法和多步法. 所谓单步法是指这类方法在计算 yn+1 时只用到前一步的 yn,这个算法的代表是 Runge-Kutta 法. Runge-Kutta 方法关于初值是稳定的, 每步的步长可以独立选取. 四阶显式 Runge-Kutta 方法是求解普通常微分方程初值问题的重要方法, 而隐式 Runge-Kutta 方法是求解刚性常微分方程初值问题的重要方法. 所谓多步法是指在计算 yn+1 时,除了用到前一步的值 yn 之外,还要用到 yn-p(p=1，2，…，k，k＞0)的值,这个算法的代表就是 Adams 方法. 1.1 龙格—库塔方法(R-K 方法) R-K 方法通过计算不同点上的函数值,并对这些函数值作线性组合,构造近似公式,利用待定系数法确定近似精度,它是非线性高阶单步法.事实上,Runge-Kutta 可以看作在(tn，tn+1)上取若干条积分曲线的若干个点的切线斜率,再进行一次(或多次)算术(或加权)平均后产生的新斜率,再按这个斜率从(tm，ym)出发,以直线带曲线向前推进一步的过 程. 由常微分方程理论,初值问题(1-1)可改写成等价的积分形式 y(t)=y0+ t t0 乙f(x，y(x))dx， 由上式及积分中值定理可得 y(tn+h)=y(tn)+hf(tn+hθ，y(tn+hθ))=yn+hk*(tn，yn，h)，0＜θ＜1. (2)但 f(tn+hθ，y(tn+hθ))未知,我们用 f 在区间［tn，tn+h］上的一些点处值的线性组合来近似它,即用下列公式代替式(2) y(tn+h)=y(tn)+h r i Σλiki, i=2，…，s, 其中 k1=fn=f(tn，yn),kj=f(tn+djh，yn+h j-1 s=1 Σβjsks)，j=2，3,…，r，dj，βjs,都是 待定的系数, 由待定系数法可以得到各阶 R-K 方法. 我们考虑 r=2 情形,应用双变量泰勒级数展式k2=f(tn+d2h，yn+β21k1h) =f(tn，yn)+h(d2 鄣 鄣t +β21k1 鄣 鄣y )f(tn，yn) +…+ h p-1 (p-1)! (d2 鄣 鄣t +β21k1 鄣 鄣y )p-1f(tn，yn)+O(hp)， 将其代入 R-K 公式,整理得 yn+1=yn+h(λ1+λ2)f(tn，yn)+h2λ2(d2ft+β21k1fy)f(tn，yn) +…+ λ2h p-1 (p-1)! (d2 鄣 鄣t +β21k1 鄣 鄣y )p-1f(tn，yn)+O(hp) 上式与 yn+1 在 yn 处的泰勒展开对比,令 h 和 h2 项的系数相同,可得 λ1+λ2=1，λ2d2= 1 2 ，λ2β21= 1 2 . d2 取不同值时可得到不同的二阶 Runge-kutta 公式. 利用上面的思想,可得经典的四阶 R-K 公式 yn+1=yn+ 1 6 h(k1+2k2+2k3+k4 ) k1=f(tn ，yn ) k2=f(tn+ 1 2 h，yn+ 1 2 hk1 ) k3=f(