matlab使用教程(31)—求解偏微分方程(2)

1.求解具有不连续性的 PDE

        此示例说明如何求解涉及物质界面的 PDE。物质界面使得问题在 x = 0 . 5 处具有不连续点，初始条件在右边界 x = 1 处具有不连续点。

        以如下分段 PDE 为例

        要在 MATLAB® 中求解该方程，您需要对方程、初始条件和边界条件编写代码，然后在调用求解器pdepe 之前选择合适的解网格。您可以将所需的函数作为局部函数包含在文件末尾（如本处所示），或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。

1.1 编写方程代码

        在编写方程代码之前，您需要确保它的形式符合 pdepe 求解器的要求。 pdepe 所需的标准形式是

现在，您可以创建一个函数以编写方程代码。该函数输入输出为
  [c,f,s] = pdex2pde(x,t,u,dudx) ：

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

function [c,f,s] = pdex2pde(x,t,u,dudx)
c = 1;
if x <= 0.5
 f = 5*dudx;
 s = -1000*exp(u);
else
 f = dudx;
 s = -exp(u);
end
end

1.2编写初始条件代码

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

初始条件为

对应的函数是

function u0 = pdex2ic(x)
if x < 1
u0 = 0;
else
u0 = 1;
end
end

1.3编写边界条件代码

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

function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)
pl = 0; 
ql = 0; 
pr = ur - 1;
qr = 0;
end

        空间网格应包括 x = 0 . 5 附近的几个值以表示不连续界面，并包括 x = 1 附近的点，因为在该点上具有不 一致的初始值 ( u 1, 0 = 1 ) 和边界值 ( u 1, t = 0 )。对于较小的 t，解的变化很快，因此请使用可以解析这 种急剧变化的时间步。

x = [0 0.1 0.2 0.3 0.4 0.45 0.475 0.5 0.525 0.55 0.6 0.7 0.8 0.9 0.95 0.975 0.99 1];
t = [0 0.001 0.005 0.01 0.05 0.1 0.5 1];

1.4 求解方程

        最后，使用对称性值 m 、PDE 方程、初始条件、边界条件以及 x 和 t 的网格来求解方程。

m = 2;
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);

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

u = sol(:,:,1);

1.5对解进行绘图

        创建在 x 和 t 的所选网格点上绘制的解 u 的曲面图。由于 m = 2 问题是在具有球面对称性的球面几何中提出的，因此解仅在径向 x 方向上变化。

surf(x,t,u)
title('Numerical solution with nonuniform mesh')
xlabel('Distance x')
ylabel('Time t')
zlabel('Solution u')

        现在，只需绘制 x 和 u 即可获得曲面图中等高线的侧视图。在 x = 0 . 5 处添加一条线，以突出材料接口的效果。

1.6 局部函数

        此处列出 PDE 求解器 pdepe 为计算解而调用的局部辅助函数。您也可以将这些函数作为它们自己的文件 保存在 MATLAB 路径上的目录中。

function [c,f,s] = pdex2pde(x,t,u,dudx) % Equation to solve
c = 1;
if x <= 0.5
 f = 5*dudx;
 s = -1000*exp(u);
else
 f = dudx;
 s = -exp(u);
end
end
%----------------------------------------------
function u0 = pdex2ic(x) %Initial conditions
if x < 1
 u0 = 0;
else
 u0 = 1;
end
end
%----------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t) % Boundary conditions
pl = 0;
ql = 0;
pr = ur - 1;
qr = 0;
end

2求解 PDE 并计算偏导数

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

        以如下 PDE 为例

        要在 MATLAB® 中求解此问题，您需要对 PDE、方程、初始条件和边界条件编写代码，然后在调用求解器 pdepe 之前选择合适的解网格。您可以将所需的函数作为局部函数包含在文件末尾（如本处所示），或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。

2.1定义物理常量

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

C.L = 1;
C.D = 0.1;
C.eta = 10;
C.K = 1;
C.Ip = 1;

2.2 编写方程代码

        在编写方程代码之前，您需要确保它的形式符合 pdepe 求解器的要求：

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

function [c,f,s] = transistorPDE(x,t,u,dudx,C)
D = C.D;
eta = C.eta;
L = C.L;
c = 1;
f = D*dudx;
s = -(D*eta/L)*dudx;
end

2.3代码初始条件

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

        初始条件为

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

2.4编写边界条件代码

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

2.5选择解网格

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

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

2.6求解方程

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

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) 处计算的解 u k 的第 k 个分量的逼近 值。对于此问题， u 只有一个分量，但通常您可以使用命令 u = sol(:,:,k) 提取第 k 个解分量。

u = sol(:,:,1);

2.7对解进行绘图

        创建在 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 即可获得曲面图中等高线的侧视图。

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

2.8 计算发射极放电电流

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

function It = serex3(t,C) % Approximate I(t) by series expansion.
Ip = C.Ip;
eta = C.eta;
D = C.D;
L = C.L;
It = 0;
for n = 1:40 % Use 40 terms
 m = (n*pi)^2 + 0.25*eta^2;
 It = It + ((n*pi)^2 / m)* exp(-(D/L^2)*m*t);
end
It = 2*Ip*((1 - exp(-eta))/eta)*It;
end

nt = length(t);
I = zeros(1,nt);
seriesI = zeros(1,nt);
iok = 2:nt;
for j = iok
 % At time t(j), compute du/dx at x = 0.
 [~,I(j)] = pdeval(m,x,u(j,:),0);
 seriesI(j) = serex3(t(j),C);
end
% Numeric solution has form I(t) = (I_p*D/K)*du(0,t)/dx
I = (C.Ip*C.D/C.K)*I;
plot(t(iok),I(iok),'o',t(iok),seriesI(iok))
legend('From PDEPE + PDEVAL','From series')
title('Emitter discharge current I(t)')
xlabel('Time t')

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

2.9局部函数

        此处列出 PDE 求解器 pdepe 为计算解而调用的局部辅助函数。您也可以将这些函数作为它们自己的文件保存在 MATLAB 路径上的目录中。

function [c,f,s] = transistorPDE(x,t,u,dudx,C) % Equation to solve
D = C.D;
eta = C.eta;
L = C.L;
c = 1;
f = D*dudx;
s = -(D*eta/L)*dudx;
end
% ----------------------------------------------------
function u0 = transistorIC(x,C) % Initial condition
K = C.K;
L = C.L;
D = C.D;
eta = C.eta;
u0 = (K*L/D)*(1 - exp(-eta*(1 - x/L)))/eta;
end
% ----------------------------------------------------
function [pl,ql,pr,qr] = transistorBC(xl,ul,xr,ur,t) % Boundary conditions
pl = ul;
ql = 0;
pr = ur;
qr = 0;
end
% ----------------------------------------------------
function It = serex3(t,C) % Approximate I(t) by series expansion.
Ip = C.Ip;
eta = C.eta;
D = C.D;
L = C.L;
It = 0;
for n = 1:40 % Use 40 terms
 m = (n*pi)^2 + 0.25*eta^2;
 It = It + ((n*pi)^2 / m)* exp(-(D/L^2)*m*t);
end
It = 2*Ip*((1 - exp(-eta))/eta)*It;
end