有限差分法
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
有限差分法
finite difference method
用差分代替微分,是有限差分法的基本出发点。是一种微分方程和积分微分方程数值解的方法。
把连续的定解区域用有限个离散点构成的网格来代替,这些离散点称作网格的节点;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似,积分用积分和来近似,于是原微分方程和定解条件就近似地代之以代数方程组,即有限差分方程组,解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。
如何根据问题的特点将定解区域作网格剖分;如何把原微分方程离散化为差分方程组以及如何解此代数方程组。此外为了保证计算过程的可行和计算结果的正确,还需从理论上分析差分方程组的性态,包括解的唯一性、存在性和差分格式的相容性、收敛性和稳定性。对于一个微分方程建立的各种差分格式,为了有实用意义,一个基本要求是它们能够任意逼近微分方程,这就是相容性要求。另外,一个差分格式是否有用,最终要看差分方程的精确解能否任意逼近微分方程的解,这就是收敛性的概念。此外,还有一个重要的概念必须考虑,即差分格式的稳定性。因为差分格式的计算过程是逐层推进的,在计算第n+1层的近似值时要用到第n层的近似值,直到与初始值有关。前面各层若有舍入误差,必然影响到后面各层的值,如果误差的影响越来越大,以致差分格式的精确解的面貌完全被掩盖,这种格式是不稳定的,相反如果误差的传播是可以控制的,就认为格式是稳定的。只有在这种情形,差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解。
最常用的方法是数值微分法,比如用差商代替微商等。另一方法叫积分插值法,因为在实际问题中得出的微分方程常常反映物理上的某种守恒原理,一般可以通过积分形式来表示。此外还可以用待定系数法构造一些精度较高的差分格式。
龙格库塔
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。对于一阶精度的欧拉公式有:
yi+1=yi+h*K1
K1=f(xi,yi)
当用点xi处的斜率近似值K1与右端点xi+1处的斜率K2的算术平均值作为平均斜率K*的近似值,那么就会得到二阶精度的改进欧拉公式:
yi+1=yi+h*( K1+ K2)/2
K1=f(xi,yi)
K2=f(xi+h,yi+h*K1)
依次类推,如果在区间[xi,xi+1]内多预估几个点上的斜率值K1、K2、……Km,并用他们的加权平均数作为平均斜率K*的近似值,显然能构造出具有很高精度的高阶计算公式。经数学推导、求解,可以得出四阶龙格-库塔公式,也就是在工程中应用广泛的经典龙格-库塔算法:
yi+1=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值。如果偏微分方程是非线性的,中心差分后得到的方程依旧是非线性方程系统,因此在时间步上推进会涉及求解非线性代数方程组。许多问题中,特别是线性扩散,代数方程中的矩阵