Mathematica 中高等微分方程的数值求解介绍

简介

MathematicaNDSolve 函数是一个通用的微分方程数值求解器. 它可以求解许多常微分方程(ODE)和一些偏微分方程(PDE). 对于常微分方程组,未知函数 的个数不限,但是,所有这些函数只能是同一个自变量 的函数. 偏微分方程可以有2个或2个以上自变量. NDSolve 也可求解一些微分-代数方程(DAE),即由微分方程和代数方程混合组成的方程组.

NDSolve[{eqn1,eqn2,...},u,{t,tmin,tmax}]求解函数 u 的数值解,t 的范围是
NDSolve[{eqn1,eqn2,...},{u1,u2,...},{t,tmin,tmax}]求解多个函数 的数值解

求解常微分方程的数值解.

NDSolveInterpolatingFunction 对象形式给出函数 的解. InterpolatingFunction 对象给出自变量 区间上, 的近似值.

通常,NDSolve 用迭代法求解. 它从特定的 值开始,然后采取一系列步骤,试着给出 整个区间的未知函数的近似值.

为了使用 NDSolve,必须给出合适的 及其导数的初始或边界条件. 这些条件给出 或许导数 在特定点 的值. 当只给出在一个点 的条件时,该方程和初始条件一起被称为初值问题. 若给出多个点 的条件时,我们称之为边值问题. NDSolve 可求解几乎所有可以写成标准形式(即对最高阶导数可解)的初值问题,但只限于线性边值问题.

以下使用 的初始条件,在 的0到2的范围内,求解 .
In[1]:=
Click for copyable input
Out[1]=

当使用 NDSolve 时,给出的初始或边界条件必须能完全确定 的解. 当使用 DSolve 求微分方程的理论解时,可指定较少的初始或边界条件. 这是因为当某个初始或边界值没有明确给出时,DSolve 会自动以任意符号常数 C[i] 来表示与初始条件相关的自由度. 但是,因为 NDSolve 必须给出数值解,所以不能表示这类附加的自由度. 因此您必须给出所需要的初始或边界条件.

一般来说,如果您有 阶导数的微分方程,那么您需要给出从0阶一直到 阶导数的初始条件或者给出 个点的边界条件.

以下是求解二阶常微分方程的初值问题,其中,在 时给出所需的两个初始条件.
In[2]:=
Click for copyable input
Out[2]=
以下是数值解的图形.
In[3]:=
Click for copyable input
Out[3]=
以下是简单的边值问题.
In[4]:=
Click for copyable input
Out[4]=

当对每个因变量求解所需的初始或边界条件都给出时,NDSolve 还可求解2个或2个以上因变量耦合在一起的微分方程组.

以下是求解一对耦合微分方程组的数值解.
In[5]:=
Click for copyable input
Out[5]=
绘出以上两个解的图形.
In[6]:=
Click for copyable input
Out[6]=

初始条件可以任何方程的形式给出,如果这些方程有多个解,NDSolve 可给出相应的多个解.

以下给出的初始条件导致多个解.
In[7]:=
Click for copyable input
Out[7]=

由于平方根函数的分支问题,NDSolve 无法给出 y'[x]=-Sqrt[y[x]^3], 的解.

以下给出 NDSolve 求出的实数部分的解(上面两个解在严格意义上的实的).
In[8]:=
Click for copyable input
Out[8]=

NDSolve 还可求解由微分方程和代数方程混合组成的微分—代数方程组(DAE). 通常,在 DAE 方程组中,有些方程并不含有导数项. 一般地说,在 DAE 方程组中导数值是不可求的,需用一种完全不同的方法来求解.

以下是一个简单的 DAE 问题.
In[9]:=
Click for copyable input
Out[9]=

请注意,尽管两方程都有导数项,但变量 没有导数项,因此 NDSolve 发出警告信息. 当用替代法将给定的微分方程组转换为一阶方程组时,其中一个方程的确变为代数方程.

另外, 因为 没有导数项,所以没有必要给出其初始条件. 事实上,给出和 DAE 方程的解相容的初始条件是相当困难的. 请参见教程 "微分代数方程的数值解" 以获得更多的信息.

以下绘出解的图形.
In[10]:=
Click for copyable input
Out[10]=

从这个图形中,可以看出 的导数很快趋于无穷大,虽然它没有明确地反映在方程组中,但这种情况意味着求解器无法继续进行.

微分方程中的未知函数并不一定要由单一的符号来表示,如果有大量的未知函数,使用 来作为函数名表示会比较方便.

以下是构建一个有 25 个耦合微分方程组和初始条件, 并对其进行求解.
In[11]:=
Click for copyable input
Out[16]=

这实际上是计算一个两端恒温的棒的热方程的近似解.(要获得更精确的解,您可以增加 .)

该解是长度为1的两端保持恒温的一维杆的热方程的近似解. 以下显示的解是作为时间函数的空间值.
In[17]:=
Click for copyable input
Out[17]=

一个未知函数也可以指定为有一个向量 (或矩阵)值. 一个未知函数的维数取决于其初始条件. 可以混合标量和矢量未知函数,只要根据 Mathematica 算法的规则,方程具有一致的维数即可. InterpolatingFunction 会给出与未知函数具有相同维数的值. 如果一个微分方程组由一个很难或不能有效用符号表示的过程决定,使用非标量变量是非常方便的.

以下使用取值为向量的未知函数来求解上面的微分方程组.
In[18]:=
Click for copyable input
Out[19]=

当指定更多的自变量时,NDSolve 还可以直接求解一些偏微分方程.

NDSolve[{eqn1,eqn2,...},u,{t,tmin,tmax},{x,xmin,xmax},...]
求解偏微分方程组,其中函数 中的 t 的范围是 x 的范围是 ,...
NDSolve[{eqn1,eqn2,...},{u1,u2,...},{t,tmin,tmax},{x,xmin,xmax},...]
求解有多个函数 的偏微分方程组

求偏微分方程组的数值解.

以下是由 NDSolve 直接得出的热方程的解.
In[20]:=
Click for copyable input
Out[20]=
以下是该解的图形.
In[21]:=
Click for copyable input
Out[21]=

NDSolve 目前使用数值直线法计算偏微分方程组的解,该方法仅限于解决至少有一个自变量的初始条件的问题. 例如,该方法不能解决椭圆偏微分方程象拉普拉斯方程,因为这些要边界值. 对于能解决的问题,直线法是相当普遍的,能处理好偏微分方程或非线性方程,而且相当快. 该方法的详细介绍请参见 "偏微分方程的数值解".

以下是数值求解具有周期边界条件的广义非线性的正弦 Gordon 方程,该方程含有2个空间变量.
In[22]:=
Click for copyable input
Out[22]=
以下是在 的结果图形.
In[23]:=
Click for copyable input
Out[23]=

如前所述,NDSolve 是通过对自变量 进行一系列步进来工作的. NDSolve 使用一个自适应的程序来确定步长. 一般来说,如果解在某个特定区域变化太快,那么 NDSolve 会减少步长以便更准确地跟踪解.

NDSolve 允许您指定所需要的精度或准确度. 一般来说,NDSolve 会采用越来越小的步长直到结果满足您所给定的 AccuracyGoal PrecisionGoal. AccuracyGoal 的设置决定了在结果中所允许的绝对误差,而PrecisionGoal 的设置决定相对误差. 如果您需要跟踪接近于零的解,那么您需要增加 AccuracyGoal 的设置. 通过设置AccuracyGoal->Infinity,你告诉 NDSolve 只使用 PrecisionGoal. 一般来说,AccuracyGoalPrecisionGoal 是用来控制在某一特定时间步骤内的误差. 对于一些微分方程,这些误差会累积,所以,最终结果的精度或准确度可能会小于您通过设置 AccuracyGoalPrecisionGoal所期望的.

NDSolve 使用您给定的 WorkingPrecision 来决定内部计算的精度. 如果您指定的 AccuracyGoalPrecisionGoal 很大,那么也需要给定一个大的 WorkingPrecision. 默认设置 Automatic 情况下, AccuracyGoalPrecisionGoalWorkingPrecision 设置的一半.

NDSolve 使用误差估计来确定所得结果是否满足规定的容许误差. 当解方程组时,它使用 NormFunction->f 选项来综合各个因变量的误差. 范数是用容许误差的大小来规范,所以 NDSolve 将尝试采取步骤使得

其中 是误差的第 个元素, 是目前解的第 个元素.

以下给出了微分方程的一个高精度解.
In[24]:=
Click for copyable input
Out[24]=
以下是解的终点值.
In[25]:=
Click for copyable input
Out[25]=

通过其自适应过程,NDSolve 能够解决"刚性"微分方程组,其中不同的因变量对于时间 有极其不同的变化率.

NDSolve 遵循一般的程序减少步长直到它准确地跟踪解. 然而,当真解有奇异点时就会有问题. 在这种情况下,NDSolve 可能会继续减少步长,永无终止. 为了避免这个问题,MaxSteps 选项指定了NDSolve会尝试求解时使用的最大步数. 对于常微分方程,缺省设置是 MaxSteps->10000.

NDSolve 在采用10,000 步后停止.
In[26]:=
Click for copyable input
Out[26]=
事实上,解在 时有一个奇异点.
In[27]:=
Click for copyable input
Out[27]=

默认设置 MaxSteps->10000 对于大多数具有平滑解的方程应该足够了. 当解是一个复杂的结构时,您可能有时不得不设置 MaxSteps 为较大的值. 当设置 MaxSteps->Infinity 时,所使用的步数将没有上限.

NDSolve 有若干不同的计算解的内置方法,同时也有增加更多的方法的机制. 默认设置 Method->Automatic 情况下,NDSolve 会选择一个适用于某微分方程组的方法. 例如,如果是刚性方程,会使用隐式法;如果方程是微分-代数方程(DAE),则会使用特殊的 DAE 方法. 在一般情况下,若没有实际求解这些问题,将无法确定微分方程解的性质:因此,缺省的 Automatic 方法可以求解很多问题,但所选择的方法对于你的问题或许并不是最好的方法. 此外,您可能想选择能保持解得某些属性的方法,如辛积分(symplectic integrators).

为一个特定系统选择一个适当的方法可能会相当困难. 让问题更为复杂的是,许多方法都有自己的设置,它们可以极大地影响解的效率和准确性. 本文档的大部分是对各种方法的描述, 以使读者对什么时候应该使用那些方法以及如何调整它们以便解决特定的问题有一个了解. 此外, NDSolve 还有一个 机制 允许您定义自己的方法,如 NDSolve 内置方法一样处理方程和结果.

NDSolve 求解时,通常有三个阶段. 首先,将给定的方程组转化为一个函数,该函数代表正规形式下方程组的右边的项; 其次,从初始条件开始以迭代方式来求解;最后,在迭代过程中储存的数据被处理成一个或多个 InterpolatingFunction 对象. 使用 中的函数您可以分别运行这些步骤,更重要的是,您可以更多地操控迭代过程. 这些步骤都是由一个 对象结合在一起,它可以保留所有解决微分方程所需要的数据.

New to Mathematica? Find your learning path »
Have a question? Ask support »