Syslab使用手册
⋮
2026a

# 求解 PDE 并计算偏导数

此示例说明如何求解一个晶体管偏微分方程 (PDE)，并使用结果获得偏导数，这是求解更大型问题的一部分。

以如下 PDE 为例

此方程出现在晶体管理论 中，u(x,t) 是描述 PNP 晶体管基极中过剩电荷载流子（或空穴）浓度的函数。D 和 η 是物理常量。该公式在区间 0≤x≤L 上对于时间 t≥0 成立。

初始条件包括常量 K，由下式给出

该问题具有由下式给出的边界条件

对于固定 x，方程 u(x,t) 的解将过剩电荷的坍塌描述为 t→∞。这种坍塌产生一种电流，称为发射极放电电流，它还有另一个常量 ：

由于在 t=0 和 t>0 时，x=0 处的边界值不一致，该公式对 t>0 有效。由于 PDE 对 u(x,t) 有闭型级数解，您可以通过解析方式和数值方式计算发射极放电电流，并对结果进行比较。

要在 Syslab 中求解此问题，您需要对 PDE、方程、初始条件和边界条件编写代码，然后在调用求解器 pdepe 之前选择合适的解网格。

# 定义物理常量

要跟踪物理常量，请创建一个NamedTuple，其中每个常量都有一个对应的字段。当您稍后为方程、初始条件和边界条件定义函数时，可以将此NamedTuple作为额外的参数传入，以便函数可以访问常量。

using TyMath
using TyPlot
C = (L=1, D=0.1, eta=10, K=1, Ip=1)

# 编写方程代码

在编写方程代码之前，您需要按照 pdepe 求解器所需的形式对其进行重写。pdepe 所需的标准形式是

以这种形式编写的 PDE 变为

因此，方程中的系数的值是

  • （没有角对称性的笛卡尔坐标）

现在，您可以创建一个函数以编写方程代码。该函数应具有签名 c,f,s = transistorPDE(x,t,u,dudx,C)：

  • x 是独立的空间变量。

  • t 是独立的时间变量。

  • u 是关于 x 和 t 微分的因变量。

  • dudx 是偏空间导数 ∂u/∂x。

  • C 是包含物理常量的额外输入。

  • 输出 c、f 和 s 对应于 pdepe 所需的标准 PDE 形式中的系数。

因此，此示例中的方程可以由以下函数表示：

function transistorPDE(x, t, u, dudx, C)
    D = C.D
    eta = C.eta
    L = C.L
    return 1, D * dudx, -(D * eta / L) * dudx
end

# 代码初始条件

接下来，编写一个返回初始条件的函数。初始条件应用在第一个时间值处，并为 x 的任何值提供 的值。使用函数签名 u0 = transistorIC(x,C) 编写函数。

初始条件为

对应的函数是

function transistorIC(x, C)
    K = C.K
    L = C.L
    D = C.D
    eta = C.eta
    return (K * L / D) * (1 - exp(-eta * (1 - x / L))) / eta
end

# 编写边界条件代码

现在，编写一个计算边界条件的函数。对于在区间 a≤x≤b 上提出的问题，边界条件应用于所有 t 以及 x=a 或 x=b。求解器所需的边界条件的标准形式是

以这种形式编写的此问题的边界条件是

对于 x=0，方程为 系数为：

同样，对于 x=1，方程为 系数为：

边界函数应使用函数 pl,ql,pr,qr = transistorBC(xl,ul,xr,ur,t)：

  • 对于左边界，输入 xl 和 ul 对应于 x 和 u。

  • 对于右边界，输入 xr 和 ur 对应于 x 和 u。

  • t 是独立的时间变量。

  • 对于左边界，输出 pl 和 ql 对应于 和 （对于此问题，x=0）。

  • 对于右边界，输出 pr 和 qr 对应于 和 （对于此问题，x=1）。

此示例中的边界条件由以下函数表示：

function transistorBC(xl, ul, xr, ur, t)
    pl = ul
    ql = 0
    pr = ur
    qr = 0
    return pl, ql, pr, qr
end

# 选择解网格

解网格定义 x 和 t 的值，求解器基于它们来计算解。由于此问题的解变化很快，请使用一个相对精细的网格，其中包含 50 个位于 0≤x≤L 区间中的空间点和 50 个位于 0≤t≤1 区间中的时间点。

x = LinRange(0, C.L, 50)
t = LinRange(0, 1, 50)

# 求解方程

最后，使用对称性值 m、PDE 方程、初始条件、边界条件以及 x 和 t 的网格来求解方程。由于 pdepe 需要 PDE 函数使用四个输入、初始条件函数使用一个输入，请创建函数，将由物理常量NamedTuple组成的作为额外输入来传入。

m = 0
eqn = (x, t, u, dudx) -> transistorPDE(x, t, u, dudx, C)
ic = (x) -> transistorIC(x, C)
sol = pdepe(m, eqn, ic, transistorBC, x, t)

pdepe 以三维数组 sol 形式返回解，其中 sol[i,j,k] 是在 t[i] 和 x[j] 处计算的解 的第 k 个分量的逼近值。sol 的大小是 length(t)×length(x)×length(u0)，因为 u0 为每个解分量指定初始条件。对于此问题，u 只有一个分量，因此 sol 是 5×20 矩阵，但通常您可以使用命令 u = sol[:,:,k] 提取第 k 个解分量。

从 sol 中提取第一个解分量。

u = sol[:,:,1];

# 对解进行绘图

创建在 x 和 t 的所选网格点上绘制的解 u 的曲面图。

surf(x,t,u)
title("Numerical Solution (50 mesh points)")
xlabel("Distance x")
ylabel("Time t")
zlabel("Solution u(x,t)")

现在，只需绘制 x 和 u 即可获得曲面图中等高线的侧视图。

figure()
plot(x,u)
xlabel("Distance x")
ylabel("Solution u(x,t)")
title("Solution profiles at several times")

# 计算发射极放电电流

使用 u(x,t) 的级数解，发射极放电电流可以表示为无穷级数:

编写一个函数，以使用级数中的 40 个项计算 I(t) 的解析解。唯一的变量是时间，但要将常量结构体指定为函数的另一个输入。

function serex3(t,C)
    Ip = C.Ip;
    eta = C.eta;
    D = C.D;
    L = C.L;

    It = 0;
    for n = 1:40
        m = (n*pi)^2 + 0.25*eta^2;
        It = It + ((n*pi)^2 / m)* exp(-(D/L^2)*m*t);
    end
    return 2*Ip*((1 - exp(-eta))/eta)*It;
end

使用 pdepe 计算的 u(x,t) 的数值解，您还可以通过以下方程计算在 x=0 处的 I(t) 的数值逼近

计算 I(t) 的解析解和数值解，并对结果绘图。使用 pdeval 计算 ∂u/∂x 在 x=0 处的值。

nt = length(t);
I = zeros(nt)
seriesI = zeros(nt)
for j in 2:nt
   _,I[j] = pdeval(m,x,u[j,:],0);
   seriesI[j] = serex3(t[j],C);
end
I = (C.Ip*C.D/C.K)*I;

plot(t[2:end],I[2:end],"o",t[2:end],seriesI[2:end])
legend("From PDEPE + PDEVAL","From series")
title("Emitter discharge current I(t)")
xlabel("Time t")

结果相当吻合。通过使用更精细的解网格，您可以进一步改进 pdepe 得出的数值结果。