第八章-常微分方程的数值解法

第八章常微分方程的数值解法,实际问题中，经常需要求解常微分方程。只有一些特殊类型的常微分方程有解析解，实际中的大多数常微分方程无法求得解析解。在多数的工程应用中，往往只需要获得若干点上的函数值，即求数值解。包含自变量、未知函数和未知函数的导数或微分的方程称为微分方程。在求解微分方程时，必须附加某种定解条件。微分方程和定解条件一起组成定解问题。定解条件通常有两种，一种是给出积分曲线在初始时刻的性态，这类条件称为初始条件，相应的定解问题就是初值问题；另一种是给出了积分曲线首末两端的性态，这类条件称为边界条件，相应的定解问题称为边值问题。,第八章常微分方程的数值解法,未知函数为一元函数的微分方程称为常微分方程，未知函数为多元函数，从而有多元函数偏导数的方程称为偏微分方程。微分方程中各阶导数的最高阶数称为微分方程的阶。本章讨论一阶常微分方程的初值问题。由常微分方程基本理论知道：当在矩形区域上连续，并且关于变量满足李普希兹（Lipschitz）条件，即对任意的，都有其中为常数，则一阶常微分方程初值问题的解存在且唯一。,第八章常微分方程的数值解法,求解常微分方程有各种各样的解析方法，但是解析方法只能求解一些特殊类型的方程，求解实际应用中得到的微分方程主要利用数值解法。对于一阶常微分方程初值问题求其数值解，就是计算出解函数在离散节点处的近似值。,第八章常微分方程的数值解法,相邻两个节点的距离称为步长，假设节点是等距的，即为常数，这时对应的函数值为这样把一个连续问题转化为一个离散型问题的过程称为离散化过程。离散化过程是把连续的微分方程初值问题转化为一个离散的差分方程问题，然后将求得的差分方程初值问题的解，作为微分方程的解在处的值的近似值。这种离散变量法又称为差分法，对应的离散方程称为差分方程或差分公式。,1欧拉（Euler）方法,欧拉公式1.数值差分法微分方程数值求解的第一步就是设法消除导数项，称为离散化。由于差分是微分的近似运算，实现离散化的基本途径是用差商代替导数。对于初值问题用向前差商代替上式中的导数项，即,1欧拉（Euler）方法,得到假设，则得到显式欧拉公式若用向后差商近似导数得到假设，则得隐式欧拉公式,1欧拉（Euler）方法,类似地，利用中心差商代替导数可得两步欧拉公式计算时，需要用到前两步信息。,1欧拉（Euler）方法,2.几何方法,1欧拉（Euler）方法,如图所示，在几何上，的精确解是一条过点的积分曲线，即图中曲线。方程的显式欧拉公式近似解对应着折线上的一系列点的纵坐标。首先过点作斜率为的直线取，得；再过点作斜率为的直线，令得。依此类推，得到折线上的一系列点，对应的纵坐标就是显式欧拉公式的数值解。,1欧拉（Euler）方法,一般地，已求得点，过此点以为斜率做直线当时，得到取。这样就可以从初始条件开始，依次计算出对应的数值解。欧拉公式的几何意义就是用一条初始点重合的折线来近似表示曲线。,1欧拉（Euler）方法,3.数值积分法把方程的两端从到积分，得到近似地计算出上式中的积分项，就能够获得的近似值；而选用不同的计算方法计算此积分项，就会得到不同的差分公式。用左矩形公式计算积分项得到就是欧拉公式。,1欧拉（Euler）方法,4.泰勒展开法对在处做二阶泰勒展开做线性近似，得到用近似值代替，则上式写为,1欧拉（Euler）方法,欧拉预估-校正方法对方程两边从到积分得利用梯形公式计算积分项有假设，则得梯形公式,1欧拉（Euler）方法,特点：隐式，不便计算；处理方法：用欧拉公式先对进行预报，记为，称为预报值，再将预报值代入梯形公式右端，得到校正的，称为欧拉预估-校正公式或改进的欧拉公式等价形式：,1欧拉（Euler）方法,例1利用欧拉公式和预估-校正公式求初值问题在区间上的数值解（取），并与精确解进行比较。解：依题意有欧拉公式：,1欧拉（Euler）方法,预估-校正公式：,1欧拉（Euler）方法,结果列表如下：,欧拉公式,预估-校正公式,精确值,1欧拉（Euler）方法,局部截断误差：假定为准确值，即，用某种数值方法计算的误差称为该数值方法计算的局部截断误差。阶方法：若某一数值方法的局部截断误差为，为正整数，则称该方法为阶方法，或说该方法具有阶精度。,1欧拉（Euler）方法,已知则欧拉方法的误差将在点泰勒展开有局部截断误差：误差为，一阶方法。,1欧拉（Euler）方法,欧拉预估-校正公式将在点泰勒展开有则,1欧拉（Euler）方法,将在点处做泰勒展开有局部截断误差：误差为，是二阶方法。小结：数值方法误差分析常用泰勒展开，展开的阶数视截断误差的阶数而定。,1欧拉（Euler）方法,例：证明下面求解公式是一个3阶方法。,1欧拉（Euler）方法,证明：,1欧拉（Euler）方法,由已知，可以分别记为局部截断误差为,1欧拉（Euler）方法,又,1欧拉（Euler）方法,1欧拉（Euler）方法,1欧拉（Euler）方法,所以公式是一个3阶方法。,1欧拉（Euler）方法,1.分别用显式欧拉公式和欧拉预估-校正公式求解初值问题在上的数值解（取）。解：欧拉公式即,1欧拉（Euler）方法,欧拉预估-校正公式即,1欧拉（Euler）方法,列表分别进行计算，结果如下,欧拉公式,欧拉预估-校正公式,2龙格-库塔（Runge-Kutta）方法,龙格-库塔方法基本思想利用微分中值定理，存在一点使得差商又，所以其中称为上的平均斜率。显然，只要提供平均斜率的一种近似，便可导出相应的计算格式。,2龙格-库塔（Runge-Kutta）方法,欧拉方法：应用点的斜率作为平均斜率的近似值改进欧拉方法：应用、两点斜率、的算术平均近似平均斜率基本思想：在区间内设法多预报几个点的斜率值，并把它们进行加权平均作为平均斜率的近似值，则可能构造精度更高的计算格式。,2龙格-库塔（Runge-Kutta）方法,二阶龙格-库塔公式在区间内任取一点，用和两点斜率的加权平均得到平均斜率，所以令其中和为待定常数。,2龙格-库塔（Runge-Kutta）方法,用欧拉公式预报的值，则则计算公式为希望通过选择和使计算格式具有较高精度。,2龙格-库塔（Runge-Kutta）方法,将在点泰勒展开有,2龙格-库塔（Runge-Kutta）方法,所以有由在点泰勒展开有比较和在的泰勒展开式可知，当时，前面公式具有二阶精度，称为二阶龙格-库塔公式。,2龙格-库塔（Runge-Kutta）方法,1）时，公式为欧拉预估-校正公式2）时，公式形式为可看作是用中点斜率作为平均斜率的近似值，称为中点公式，具有二阶精度。,2龙格-库塔（Runge-Kutta）方法,三阶龙格-库塔公式截断误差为。,2龙格-库塔（Runge-Kutta）方法,四阶龙格-库塔公式截断误差为。,2龙格-库塔（Runge-Kutta）方法,举例：随机共振现象根据随机共振理论，当淹没在强背景噪声中的微弱周期信号通过某些特殊非线性系统时，在噪声与信号的协同作用下，会发生噪声能量向信号能量的转移，使得非线性系统的输出信噪比得到增强。随机共振系统包括三个必不可少的因素：非线性系统、输入信号和噪声。最简单的非线性系统为郎之万方程描述的双稳态系统，最简单的信号和噪声为单频的周期信号和高斯白噪声。,2龙格-库塔（Runge-Kutta）方法,双稳态系统模型为式中为外加周期信号，是均值为0、强度为D的白噪声，满足白噪声可以表示为式中为噪声的方差，是均值为0、方差为1的白噪声。记,2龙格-库塔（Runge-Kutta）方法,采用四阶龙格-库塔公式进行求解，公式为,2龙格-库塔（Runge-Kutta）方法,t=0:40000;n=randn(1,40001)*0.04;%强度为0.0008的高斯白噪声fs=10000;%采样频率s=0.003*cos(2*pi*t/fs);%输入周期信号p=s+n;%输入信号+噪声pp=fft(p);figure(1);plot(s);gridon;title(输入周期信号);figure(2);plot(p);gridon;title(输入周期信号+噪声);figure(3);plot(abs(pp);gridon;title(系统输入频谱);a=0.01;b=0.01;h=1;x=zeros(1,40001);fori=1:40000k1=h*(a*x(i)-b*(x(i)3)+p(i);k2=h*(a*(x(i)+k1/2)-b*(x(i)+k1/2)3)+p(i);k3=h*(a*(x(i)+k2/2)-b*(x(i)+k2/2)3)+p(i+1);k4=h*(a*(x(i)+k3)-b*(x(i)+k3)3)+p(i+1);x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6;Endy=fft(x);figure(4);plot(x);gridon;title(SR处理后信号);figure(5);plot(abs(y);gridon;title(SR处理后信号频谱);,2龙格-库塔（Runge-Kutta