偏微分方程数值详讲

3年前浏览3028

过冷水诚挚邀请你加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载，qq：927550334

一般而言,偏微分方程解析求解只有在定解问题比较简单,求解区域比较规则条件下才有可能应用，可解决问题的局限性比较大。实际工程问题中遇到的情况都比较复杂,或由于边界不规则,解析解方法比较困难,大多数问题需采用数值解法。偏微分方程有多种数值解法,如有限差分法、正交配置法、线上法(MOL)、有限元法等，本期过冷水就和大家详细了解一下数值方法的具体原理和实用案例。

线上法指将偏微分方程仅有一个自变量保持连续,而将其它变量进行离散，从而把偏微分方程转变为常微分方程组进行求解的数值方法。通常线上法选择离散的自变量是空间变量,而将时间变量保持连续,离散方法可以选择有限差分法或配置法。

考虑一维扩散方程:

其中扩散系数D为常数。将空间变量x采用差分法进行离散。令:

则：

由于转化过程中,原方程的边界条件被包含在离散化过程中,原方程变为常微分方程组的初值问题求解。

对于二维过程的离散也类似,例如,二维扩散方程;

离散：

讲完数值解的离散方法后，求解偏微分方程的具体实现程序

Pdepe函数求解偏微分方程

pdepe函数用于求解如下类型问题：

满足条件：

pdepe函数调用格式

sol= pdepe(m，pdepe，icfun，bcfun，xmesh， tspan)
sol = pdepe(m，pdepe，icfun，bcfun，xmesh，tspan，options)
sol = pdepe(m，pdepe，icfun，bcfun，xmesh，tspan，options，p1，p2，...)

( 1) m就是式(7,10)级数, m=0，1，2分别代表平板、圆柱和球形。

(2) pdefun是描述偏微分方程的函数,其格式为[c，f，s]=pdefun( x，t，u，dudx)。输出变量c，f,s是列向量,输人变量x和 t为标量,表示空间和时间的自变量,u和dudx分别是问题的解u(x，t)和它对x的偏导数的近似。

(3) icfun描述定解问题初始条件的函数,其格式为u=icfun(x)。icfun计算和返回解的初始值。

( 4) bcfun是描述定解问题边界条件的函数,格式为[pl，ql，pr，qr]= bcfun( xl，ul,xr，ur，t)。其中ul是在左边界xl=α处的近似解,ur是在右边界xr=b处的近似解。pl和ql是p和q在xl处的列向量值,同样pr和qr是p和q在xr处的列向量值。当m>0和 a=0时﹐解的有界性要求通量f在a=0处为0, pdepe 会自动处理,并忽略pl和ql的值。

(5) xmesh是空间变量x的网点向量, pdepe不会自动选择。xmesh=[x0，xl，…,xn],满足 x0<x1<…<xn,且 x....的长度必须大于3, xmesh(1)=a和 xmesh(end)=b.。pdepe的求解效率与x..的选择好坏关系很大。通常在梯度较大处应加密网格。

(6) tspan是时间变量t 的网点向量,是用户在该处想得到数值解的时间向量。tspan—[t0，t1，…，tn],满足t0<tl<…<tf,且 t.span的长度必须大于3，tspan( 1)=t0和tspan(end)=tf。求解效率和 tspan的疏密关系不大。

(7) sol是一个三维数组,存放数值解。sol(:，:，i)=ui是解向量的第i个分量。sol(j，k，i)=u(j，k)是解向量ui在(x，t)=(tspan(j)，xmesh(k))处的数值解。sol(j，:.i)是解向量第i个分向量ui在时间 t.span(j)和网格点xmesh(:)处的数值解。

(8) options是 MATLAB常微分方程求解函数( pdepe函数需调用ode15s进行时间积分)options 的一部分,如RelTol，AbsTol，NormControl，InitiaIStep 和 MaxStep。一般不用调整options 的值即可获得满意的解,如果需要改变则可以采用odeset 函数,这与常微分方程求解是相同的。

(9) p1.p2，...是通过pdepe传递给pdefun，icfun和 bcfun的参数。

案例演示

对于一维扩散问题

边界条件：

初始条件：

当D=1时,试求解以上方程,采用图形将计算结果输出,并与x = 0.3，t = 0.005、0.01、0.02的精确解0.5966、0.5799、0.5334比较。

根据函数可求得图像为：

pdepe求解程序为:

function  pde
m=0;
xmesh= 0:0.01:1;
tspan=0:0.001:0.1;
sol = pdepe(m,@pdefun,@icfun,@bcfun,xmesh,tspan);
surf( xmesh,tspan,sol)
xlabel('x'), ylabel('time'), zlabel('w');
shading interp
figure
plot( tspan,sol(1:20:101,:))
xlabel('time')
ylabel('w')
legend('x=0','x=0.2','x=0.4','x=0.6','x=0.8','x=1.0')
disp('The anlytic results for t=0.005 0.01 0.02 @x=0.3 are;' )
disp([0.5966 0.5799 0.5334])
disp('The corresponding results obtained from pdepe are:')
disp([sol(6,31),sol( 11,31),sol(21,31)])
function [c,f,s]= pdefun(x,t,u,du)
    c= 1;
    D=1;
    f=D*du;
    s=0;
end
function u=icfun(x)
    if  x>=0&x<=1/2
        u=2*x;
    elseif x>=1/2&x<=1
        u=2*(1-x);
    end
end
function [pl,ql,pr,qr]= bcfun(xl,ul,xr,ur,t)
    pl = ul;
    ql =0;
    pr = ur;
    qr=0;
end
end

该结果相对来说还算是满意的，在调试过程中过冷水发现在偏微分计算中函数调用形式和以前接触的方法有点不同，使用起来不太习惯，调试花费了好久时间，看来知识还需要精进.