MATLAB求解偏微分方程（PDE）教程
MATLAB求解偏微分方程（PDE）教程
本文是一篇关于使用MATLAB求解偏微分方程（PDE）的教程，详细介绍了两个具体的案例：求解具有不连续性的PDE和求解晶体管PDE并计算偏导数。每个案例都包含了完整的代码示例，包括方程、初始条件、边界条件的编写，以及如何使用MATLAB的pdepe求解器进行求解。
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,t0 )的 值。使用函数签名**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 附近的点，因为在该点上具有不 一致的初始值 (u1, 0 = 1) 和边界值 (u1,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)处计算的解uk的第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,t0） 的值。使用函数签名**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)处计算的解uk的第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