常微分方程数值解.ppt

Euler方法的局部截断误差： 故Tn+1= O(h2)，p=1， （设yn=y(xn)） 其中 称局部截断误差主项。 即Euler方法具1阶精度。 （设yn=y(xn)） 故Tn+1= O(h3)，p=2， 梯形方法的局部截断误差： 局部截断误差主项为： 梯形方法具2阶精度。 6.3 Runge－Kutta方法 一、Runge-Kutta方法的基本思想 由Taylor展式 Tn+1= O(hp+1)，若提高p，可提高精度。 但因 …… 高阶导数计算复杂，故可从另外角度考虑。 分析Euler公式及改进的Euler公式： 局部截断误差：O(h2) 局部截断误差：O(h3) 可用f(x,y)在某些点处值的线性组合得yn+1，增加计算f(x,y)的次数可提高阶数。 设法计算f(x,y)在某些点上的函数值，然后对这些函数值作线性组合，构造近似计算公式，再把近似公式和解的泰勒展开式相比较，使前面的若干项吻合，从而获得达到一定精度的数值计算公式 。 Runge-Kutta方法的基本思想： 设 ci , λi , μij 为待定常数。 上面第一个式子的右端在(xn,yn)作泰勒展开后,按h的幂次作升序排列 ： 再与初值问题的精确解y(x)在点x=xn处的泰勒展开式 相比较，使其有尽可能多的项重合。 例如，要求 就得到p个方程,从而定出参数ci ，?i，?ij ，再代入K1，K2，…， Kr的表达式，就可得到计算微分方程初值问题的数值计算公式 ： 若 Tn+1= O(hp+1)，则称其为p 阶r 级R－K方法。 上式称为r 级Runge－Kutta方法的计算公式。 当r=1时，就是Euler方法。 要使Runge-Kutta公式具有更高的阶p，就要增加r 的值。下面我们只就r =2推导R-K方法。 二、二阶Runge-Kutta方法 其中 c1, c2, ?2, ?21 待定。 上式的局部截断误差为： 又 由 利用二元函数的Taylor展开，得 代入Tn+1的表达式，得 即 c1 = 1-a ， ?2 = ?21 =1/(2a) 要使上式p=2阶，则需 方程组解不唯一，可令c2=a ? 0 ，则 满足上述条件的公式都为2阶R-K公式。 称中点公式，相当于数值积分的中矩形公式： 如取a= ? ，则c1= c2= ?， ?2=?21 =1，即为改进Euler公式。 若取a= 1，则c1= 0，c2= 1， ?2=?21 = ? ，得 例:蛇形曲线的初值问题 令f(x,y)=y/x –2y2, 取f(0,0)=1, h=0.2, xn=nh , ( n = 1,2,…,15) 2阶龙格-库塔公式计算格式： k1=yn/xn – 2yn2, k2 = (yn+hk1)/(xn+h) – 2(yn+hk1)2 yn+1=yn + 0.5h[ k1 + k2] x0=0;y0=0;h=.2; x=.2:h:3; k1=1; k2=(y0+h*k1)/x(1)-2*(y0+h*k1)^2; y(1)=y0+.5*h*(k1+k2); for n=1:14 k1=y(n)/x(n)-2*y(n)^2; k2=(y(n)+h*k1)/x(n+1)-2*(y(n)+h*k1)^2; y(n+1)=y(n)+0.5*h*(k1+k2); end y1=x./(1+x.^2); subplot(221) plot(x,y1,x,y1,b*) subplot(222) plot(x,y,x,y,o) subplot(223) plot(x,y,x,y,o,x,y1,x,y1,b*‘) 三、三阶与四阶Runge-Kutta方法 当r=3时，R-K公式表示为 其中 为8个待定常数。 上式的局部截断误差为 类似二阶的推导过程，将K2, K3按二元函数展开，使Tn+1= O(h4)，得 方程有8个未知数，解不唯一。 满足该条件的公式统称为三阶R-K公式。 其中一个常用公式为： 当r=4时，利用相同的推导过程，经过较复杂的计算，可以得出四阶R-K公式的成立条件。 下列经典公式是其中常用的一个： 四、一阶常微分方程组和高阶微分方程的 数值解法简介 一阶常微分方程组的数值解法： 下列包含多个一阶常微分方程的初值问题： 称为一阶常微分方程组的初值问题。 引进向量记号: 则上述一阶常微分方程组的初值问题化为矩阵形式： 它在形式上跟单个微分方程的初值问题形式完全相同，只是函数变成了向量函数。故前面介绍的一切数值方法都适用，只要把函数换成向量函数即可。 k1=f(xn,yn), k2=f(xn+0.5h,yn+0.5hk1) k3=f(xn+0.5h,yn+0