二阶pde数值求解算法

Orbifold

Last updated on Jan 24, 2024 4 min read

二阶偏微分方程分类

$A U_{x x} + 2 B U_{x y} + C U_{y y} + D U_{x} + E U_{y} + F U = 0$

若 $B^{2} - A C < 0$ ,椭圆型；
若 $B^{2} - A C = 0$ ,抛物线型；
若 $B^{2} - A C > 0$ ,双曲型；

边界条件

第一类边界条件，给定边界 $Γ$ 上的函数值 $ϕ |_{Γ}$ ;
第二类边界条件，给定边界 $Γ$ 上的函数导数值 $\frac{\partial ϕ}{\partial n} |_{Γ}$ ;
第三类边界条件，给定边界 $Γ$ 上的函数值与导数值线性组合 $a ϕ + b \frac{\partial ϕ}{\partial n} |_{Γ}$ ;

椭圆形偏微分方程

主要考虑Poisson方程和Laplace方程:

$- [\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}] ϕ = S (x, y)$

区域离散化：

$- [ϕ_{i + 1, j} + ϕ_{i - 1, j} + ϕ_{i, j + 1} + ϕ_{i, j - 1} - 4 ϕ_{i j}] = h^{2} S_{i j}$

或表示为：

$- [(δ_{i}^{2} ϕ)_{i j} + (δ_{j}^{2} ϕ)_{i j}] = h^{2} S_{i j}$

$\begin{array}{l} (δ_{i}^{2} ϕ)_{i j} \equiv ϕ_{i + 1, j} + ϕ_{i - 1, j} - 2 ϕ_{i j} (δ_{j}^{2} ϕ)_{i j} \equiv ϕ_{i, j + 1} + ϕ_{i, j - 1} - 2 ϕ_{i j} \end{array}$

正方形边界代码示例:

% 边值条件与初值选取
xa=0;xb=pi;
ya=0;yb=pi;
N=100;
M=100;
h1=(xb-xa)/N;
h2=(yb-ya)/M;
xinit=xa:h1:xb;
yinit=ya:h2:yb;
u=zeros(N+1,M+1);
u(:,M+1)=sin(xinit);
uold=u;
u(2:N,2:M)=1/(2*pi);% 边界平均值为1/(2*pi)以初始化内点

Jacobi迭代法

$φ_{i j}^{n + 1} = \frac{1}{4} [φ_{i - 1, j}^{n} + φ_{i + 1, j}^{n} + φ_{i, j - 1}^{n} + φ_{i, j + 1}^{n} + h^{2} S_{i j}^{n}]$

% 正方形边界问题下的jacobi迭代法实现
function usolve=jacobi_iteration(u,s,h,uold,N,M,err)
if nargin <5         
    err = 10^-6;
end  
iter_num=0;
while max(max(abs(u-uold)))>=err
      uold=u;
      u(2:N,2:M)=0.25*(u(3:N+1,2:M)+u(1:N-1,2:M)+u(2:N,3:M+1)+u(2:N,1:M-1)+h^2*s(2:N,2:M));
      iter_num=iter_num+1;
end
usolve=u;
% fprintf("The error of Jacobi method is %f\n",max(max(abs(u-uold))))
fprintf("The iteration number of Jacobi method is：%d\n",iter_num)
end

Gauss-Seidal法

$φ_{i j}^{n + 1} = (1 - ω) φ_{i j}^{n} + \frac{ω}{4} [φ_{i - 1, j}^{n + 1} + φ_{i + 1, j}^{n} + φ_{i, j - 1}^{n + 1} + φ_{i, j + 1}^{n} + h^{2} S_{i j}^{n}]$

其中 $0 < ω < 2$ ，最佳值需要依据情况确定。

% 正方形边界问题下的Gauss_Seidel迭代法实现
function usolve=Gauss_Seidel(u,S,h,w,N,M,err)
if nargin <5         
    err = 10^-8;
end    
iter_num=0;
merr=err+1;
while merr>err
    iter_num= iter_num+1;
    merr=0;
        for i=2:N
            for j=2:M
                %Gauss-Seidel迭代法
                relax=w/4*(u(i,j+1)+u(i,j-1)+u(i+1,j)+u(i-1,j)-4*u(i,j)+h^2*S(i,j));
                u(i,j)=u(i,j)+relax;
                if abs(relax)>merr
                    merr=abs(relax);
                end
            end
        end
end
 usolve=u;       
%  fprintf("The error of Gauss-Seidel method is%f\n",abs(merr))
 fprintf("The iteration number of Gauss-Seidel method is：%d\n",iter_num)
end

抛物线形偏微分方程

主要以扩散方程和含时薛定谔方程为主：

一维扩散方程：

$\frac{\partial Φ}{\partial t} = \frac{\partial^{2} Φ}{\partial x^{2}} + S (x, t)$

一维扩散问题边界示例

%格点化条件
L=1;
t0=1;
a=1;%u_t=a^2u_xx
dx=0.05;x=0:dx:L;N=length(x);
dt=0.001;t=0:dt:t0;M=length(t);
[X,T]=meshgrid(x,t);
y=zeros(N,M);
y(1:(N+1)/2,1) = 2*x(1:(N+1)/2); % 边界条件
y((N+1)/2:N,1) = 2*(1-x((N+1)/2:N));

显式差分

$Φ_{i}^{n + 1} = r Φ_{i + 1}^{n} + (1 - 2 r) Φ_{i}^{n} + r Φ_{i - 1}^{n} + S_{i}^{n} Δ t, r = \frac{Δ t}{h^{2}}$

只有当 $r < \frac{1}{2}$ 时算法才稳定。

若记

(δ^{2} Φ^{n})_{i} = Φ_{i + 1}^{n} - 2 Φ_{i}^{n} + Φ_{i - 1}^{n}

则显式差分也可写成：

$Φ^{n + 1} = (1 - H Δ t) Φ^{n} + S^{n} Δ t$

其中：

$H Φ_{i} \equiv - \frac{1}{h^{2}} (δ^{2} Φ)_{i}$

%----显式差分-------
function  uforward=forward(y,dx,dt,M,N,a)
rr=dt/dx^2*a;
for j=1:M-1    %时间上的步数
    for i=2:N-1   %空间上的点
        y(i,j+1)=rr*y(i+1,j)+(1-2*rr)*y(i,j)+rr*y(i-1,j);
    end
end
uforward=y';
end

直接隐式差分

[\begin{array}{l} 1 + 2 r & - r \\ - r & 1 + 2 r & - r \\ ⋱ & ⋱ & ⋱ \\ - r & 1 + 2 r & - r \\ - r & 1 + 2 r \end{array}] [\begin{array}{l} Φ_{1}^{n + 1} \\ Φ_{2}^{n + 1} \\ \dots \\ Φ_{N - 2}^{n + 1} \\ Φ_{N - 1}^{n + 1} \end{array}] = [\begin{array}{l} Φ_{1}^{n} + S_{1}^{n} Δ t + r Φ_{0}^{n + 1} \\ Φ_{2}^{n} + S_{2}^{n} Δ t \\ . . . \\ Φ_{N - 2}^{n} + S_{N - 2}^{n} Δ t \\ Φ_{N - 1}^{n} + S_{N - 1}^{n} Δ t + r Φ_{N}^{n + 1} \end{array}]

用算符表示为：

$Φ^{n + 1} = \frac{1}{1 + H Δ t} [Φ^{n} + S^{n} Δ t]$

%----直接隐式差分----
function  uback=back(y,dx,dt,M,N,a)
rr=dt/dx^2*a;
A=diag((1+2*rr)*ones(N-2,1),0) + diag(-rr*ones(N-3,1),-1) + diag(-rr*ones(N-3,1),1);
for j=1:M-1    %时间上的步数
    for i=2:N-1   %空间上的点
        b = y(2:N-1,j);
        b(1) = b(1) + rr*y(1,j+1);
        b(N-2) = b(N-2) + rr*y(N,j+1);
        y(2:N-1,j+1) = soldiag(A,b);
    end
end
uback=y';
end
%--soldiag---
function x = soldiag(A,d)
%求解Ax=d，A为系数矩阵
    n = length(d);
    x = zeros(n,1);
    for i = 2:n
        w = A(i,i-1)/A(i-1,i-1);
        A(i,i) = A(i,i)-w*A(i-1,i);
        d(i) = d(i)-w*d(i-1);
    end
    x(n) = d(n)/A(n,n);
    for i = (n-1):(-1):1
        x(i) = (d(i)-A(i,i+1)*x(i+1))/A(i,i);
    end
end

平均隐式差分

[\begin{array}{l} 1 + r & - r / 2 \\ - r / 2 & 1 + r & - r / 2 \\ ⋱ & ⋱ & ⋱ \\ - r / 2 & 1 + r & - r / 2 \\ - r / 2 & 1 + r \end{array}] [\begin{array}{l} Φ_{1}^{n + 1} \\ Φ_{2}^{n + 1} \\ \dots \\ Φ_{N - 2}^{n + 1} \\ Φ_{N - 1}^{n + 1} \end{array}] = [\begin{array}{l} \frac{r}{2} Φ_{2}^{n} + (1 - r) Φ_{1}^{n} + \frac{r}{2} Φ_{0}^{n} + S_{0}^{n} Δ t + \frac{r}{2} Φ_{0}^{n + 1} \\ \frac{r}{2} Φ_{3}^{n} + (1 - r) Φ_{2}^{n} + \frac{r}{2} Φ_{1}^{n} + S_{1}^{n} Δ t \\ . . . \\ \frac{r}{2} Φ_{N - 1}^{n} + (1 - r) Φ_{N - 2}^{n} + \frac{r}{2} Φ_{N - 3}^{n} + S_{N - 2}^{n} Δ t \\ \frac{r}{2} Φ_{N}^{n} + (1 - r) Φ_{N - 1}^{n} + \frac{r}{2} Φ_{N - 2}^{n} + S_{N - 1}^{n} Δ t + \frac{r}{2} Φ_{N}^{n + 1} \end{array}]

用算符表示为：

$Φ^{n + 1} = \frac{1}{1 + \frac{1}{2} H Δ t} [(1 - \frac{1}{2} H Δ t) Φ^{n} + S^{n} Δ t]$

%---------平均隐式差分-----
function uck=Crank_Nicolson(y,dx,dt,M,N,a)
rr=dt/dx^2*a;
A = diag((1+rr)*ones(N-2,1),0) + diag(-rr/2*ones(N-3,1),-1) + diag(-rr/2*ones(N-3,1),1);
for j = 1:M-1
    b = (1-rr)*y(2:N-1,j) + rr/2*y(1:N-2,j) + rr/2*y(3:N,j);
    b(1) = b(1) + rr/2*y(1,j+1);
    b(N-2) = b(N-2) + rr/2*y(N,j+1);
    y(2:N-1,j+1) = soldiag(A,b);
end
uck = y';
end

一维薛定谔方程的平均隐式差分

$i ℏ \frac{\partial ϕ}{\partial t} = - \frac{ℏ^{2}}{2 m} \nabla^{2} ϕ + V ϕ$

平均隐式差分为：

$ϕ^{n + 1} = (\frac{1 - i \frac{1}{2} H Δ t}{1 + i \frac{1}{2} H Δ t}) ϕ^{n}, H = - \frac{\partial^{2}}{\partial x^{2}} + V$

实际计算时：

$ϕ_{}^{n + 1} = (\frac{2}{1 + i \frac{1}{2} H Δ t} - 1) ϕ_{}^{n} \equiv χ - ϕ_{}^{n}$

其中 $χ$ 满足：

$χ_{j - 1} + [- 2 + \frac{2 i h^{2}}{Δ t} - h^{2} V_{j}] χ_{j} + χ_{j + 1} = \frac{4 i h^{2}}{Δ t} ϕ_{j}^{n}, j = 1, \cdot \cdot \cdot, N - 1$

或矩阵形式：

令：

$r = - 2 + \frac{2 i h^{2}}{Δ t} - h^{2} V_{j}$

[\begin{array}{l} r & 1 \\ 1 & r & 1 \\ ⋱ & ⋱ & ⋱ \\ 1 & r & 1 \\ 1 & r \end{array}] [\begin{array}{l} χ_{1} \\ χ_{2} \\ \dots \\ χ_{N - 2} \\ χ_{N - 1} \end{array}] = [\begin{array}{l} \frac{4 i h^{2}}{Δ t} ϕ_{1}^{n} - χ_{0} \\ \frac{4 i h^{2}}{Δ t} ϕ_{2}^{n} \\ . . . \\ \frac{4 i h^{2}}{Δ t} ϕ_{N - 2}^{n} \\ \frac{4 i h^{2}}{Δ t} ϕ_{N - 1}^{n} + χ_{N} \end{array}]

一维高斯波包初值边界值示例

%设定条件
xa=-30;xb=30;dx=0.1;
x=xa:dx:xb;
ta=0;tb=5;dt=0.01;t=ta:dt:tb;
N = length(x);
M = length(t);
V = transpose(arrayfun(@v,x));
sigma=3;k0=2.5;%Gauss波包参数
x0=-10;              %Gauss波包参数

psi0 = exp(k0*1i.*x).*exp(-(x-x0).^2*log10(2)/(2*sigma^2)); 
psi0(1) = 0; psi0(N) = 0; % 波函数边界条件
psi0 = transpose(psi0);

平均隐式差分

%平均隐式差分
psi=zeros(N,M);
psi(:,1)=psi0;
chi = zeros(N,M-1);
r = 1i*dt/(dx^2);
A = diag((1+r)*ones(N-2,1)+1i/2*dt*V(2:N-1)) + diag(-r/2*ones(N-3,1),-1) + diag(-r/2*ones(N-3,1),1);
for j = 1:M-1
    b = 2*psi(2:N-1,j);
    b(1) = b(1) + r/2*chi(1,j);
    b(N-2) = b(N-2) + r/2*chi(N,j);
    chi(2:N-1,j) = inv(A)b;
    psi(2:N-1,j+1) = chi(2:N-1,j)-psi(2:N-1,j);
end

计算物理

二阶pde数值求解算法

二阶偏微分方程分类

边界条件

椭圆形偏微分方程

Jacobi迭代法

Gauss-Seidal法

抛物线形偏微分方程

显式差分

直接隐式差分

平均隐式差分

一维薛定谔方程的平均隐式差分

Orbifold

Undergraduated Student