有限差分法

有限差分法
finite difference method
用差分代替微分，是有限差分法的基本出发点。

是一种微分方程和积分微分方程数值解的方法。

把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。

然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。

如何根据问题的特点将定解区域作网格剖分；如何把原微分方程离散化为差分方程组以及如何解此代数方程组。

此外为了保证计算过程的可行和计算结果的正确，还需从理论上分析差分方程组的性态，包括解的唯一性、存在性和差分格式的相容性、收敛性和稳定性。

对于一个微分方程建立的各种差分格式，为了有实用意义，一个基本要求是它们能够任意逼近微分方程，这就是相容性要求。

另外，一个差分格式是否有用，最终要看差分方程的精确解能否任意逼近微分方程的解，这就是收敛性的概念。

此外，还有一个重要的概念必须考虑，即差分格式的稳定性。

因为差分格式的计算过程是逐层推进的，在计算第n＋1层的近似值时要用到第n层的近似值，直到与初始值有关。

前面各层若有舍入误差，必然影响到后面各层的值，如果误差的影响越来越大，以致差分格式的精确解的面貌完全被掩盖，这种格式是不稳定的，相反如果误差的传播是可以控制的，就认为格式是稳定的。

只有在这种情形，差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解。

最常用的方法是数值微分法，比如用差商代替微商等。

另一方法叫积分插值法，因为在实际问题中得出的微分方程常常反映物理上的某种守恒原理，一般可以通过积分形式来表示。

此外还可以用待定系数法构造一些精度较高的差分格式。

龙格库塔
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。

由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。

该算法是构建在数学支持的基础之上的。

对于一阶精度的欧拉公式有：
yi＋１＝yi＋h*K1
K1=f(xi,yi)
当用点xi处的斜率近似值K1与右端点xi＋１处的斜率K2的算术平均值作为平均斜率K*的近似值，那么就会得到二阶精度的改进欧拉公式：
yi＋１＝yi＋h*( K1+ K2)/2
K1=f(xi,yi)
K2=f(xi+h,yi+h*K1)
依次类推，如果在区间[xi,xi+1]内多预估几个点上的斜率值K1、K2、……Km，并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。

经数学推导、求解，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法：
yi＋１＝yi＋h*( K1+ 2*K2 +2*K3+ K4)/6
K1=f(xi,yi)
K2=f(xi+h/2,yi+h*K1/2)
K3=f(xi+h/2,yi+h*K2/2)
K4=f(xi+h,yi+h*K3)
ode45:使用四阶龙格-库塔法求微分方程，调用格式为：
[t，y]=ode45(‘fname’,tspan,y0)
其中fname为由M函数定义的线性或者非线性微分方程的句柄函数名。

tspan形式为[t0,tf]，表示求解区间。

y0是初始状态列向量。

t和y分别给出时间向量和相应的状态向量。

四阶龙格-库塔法(ode45)：计算原理为预报-校正法，预报值采用Euler算出，下式又作了3次校正，因此计算精度远高于Euler算法。

ode23:使用二阶龙格-库塔法求微分方程，调用格式为：
[t，y]=ode23(‘fname’,tspan,y0)
eg. 二阶非线性系统的微分方程：
x″ + 0.5*x′+ 2*x + x^2 = 0
求系统在初始条件为x(0)=1，x′(0)=0的数值解。

建立M函数：
function xdot=odetest(t,x)
%龙格-库塔算法测试
%二阶非线性系统微分方程xdotdot + 0.5*xdot + 2*x + x^2 = 0
%求系统在初始条件为x(0)=1,xdot(0)=0的数值解
xdot=zeros(2,1);
xdot(1)=x(2);
xdot(2)=-0.5*x(2)-2*x(1)-x(1)^2;
tspam=[0 20];
x0=[0 1];
[t,x]=ode23(@odetest,tspan,x0);
plot(t,x)
微分方程数值解曲线如下图：
Crank-Nicolson 方法
数值分析中，Crank-Nicolson方法是有限差分方法中的一种，用于数值求解热方程以及形式类似的偏微分方程。

它在时间方向上是隐式的二阶方法，数值稳定。

该方法诞生于20世纪，由John Crank与Phyllis Nicolson发展。

对于扩散方程（包括许多其他方程），可以证明Crank-Nicolson方法无条件稳定。

但是，如果时间步长与空间步长平方的比值过大（一般地，大于1/2），近似解中将存在虚假的振荡或衰减。

基于这个原因，当要求大时间步或高空间分辨率的时候，往往会采用数值精确较差的后向欧拉方法进行计算，这样即可以保证稳定，又避免了解的伪振荡。

目录
[隐藏]
•1方法
•2示例
o 2.1线性扩散问题
o 2.2一维多通道连接的扩散问题
[编辑]方法
Crank-Nicolson方法在空间域上的使用中心差分；而时间域上应用梯形公式，保证了时间域上的二阶收敛。

例如，一维偏微分方程
令，则通过Crank-Nicolson方法导出的差分方程是第n步上采用前向欧拉方法与第n+1步上采用后向欧拉方法的平均值（注意，Crank-Nicolson方法本身不是这两种方法简单地取平均，方程对解隐式依赖）。

(前向欧拉方法)
(后向欧拉方法)
(Crank-Nicolson)
对于F，通过中心差分方法使其在空间上是离散的。

注意，这是一个隐式方法，需要求解代数方程组以得到时间域上的下一个u值。

如果偏微分方程是非线性的，中心差分后得到的方程依旧是非线性方程系统，因此在时间步上推进会涉及求解非线性代数方程组。

许多问题中，特别是线性扩散，代数方程中的矩阵
是三对角的，通过三对角矩阵算法可以高效求解，这样，算法的时间复杂度由直接求解全矩阵的转化为。

[编辑]示例
[编辑]线性扩散问题
•线性扩散方程
通过Crank-Nicolson方法将得到离散方程
引入变量:
这是一个三对角问题，应用三对角矩阵算法（追赶法）即可得到，而不需要对矩阵直接求逆。

•准线性扩散方程
离散化后则会得到非线性方程系统。

但是某些情况下，通过使用a的旧值，即用替代，可将问题线性化。

其他时候，也可能在保证稳定性的基础上使用显式方法估计
[编辑]一维多通道连接的扩散问题
这种模型可以用于描述水流中含稳定污染流，但只有一维信息的情况。

它可以简化为一维问题并得到有价值的信息。

我们对水中污染溶质富集的问题进行建模，这种问题由三部分组成：已知的扩散方程（D x为常量），平流分量（即由速度场导致的系统在空间上的变化，表示为常量Ux），以及与纵向通道k旁流的相互作用。

其中C表示污染物的富集水平，下标N和M分别对应上一通道和下一通道。

Crank-Nicolson方法（i对应位置，j对应时间）将以上偏微分方程中的每个部分变换为现在引入以下常量用于简化计算：
把<1>, <2>, <3>, <4>, <5>, <6>, α, β和λ代入<0>. 把新时间项(j+1)代入到左边，当前时间项(j)代入到右边，将得到
第一个通道只能与下一个通道(M)有关系，因此表达式可以简化为：
同样地，最后一个通道只与前一个通道（N）有关联，因此表达式可以简化为
为求解此线性方程组，我们需要知道边界条件在通道始端就已经给定了。

: 当前时间步某通道的初始条件
: 下一时间步某通道的初始条件
: 前一通道到当前时间步下某通道的初始条件
: 下一通道到当前时间步下某通道的初始条件
对于通道的末端最后一个节点，最方便的条件是是绝热近似，则
当且只当
时，这一条件才被满足。

以3个通道，5个节点为例，可以将线性系统问题表示为
其中，
需要清楚的是，AA和BB是由四个不同子矩阵组成的矩阵，
其中上述矩阵的的矩阵元对应于下一个矩阵和额外的4x4零矩阵。

请注意，矩阵AA和BB的大小为12x12
&。