Matlab的偏微分方程工具箱求解方法

这一节我们主要用matlab自带的偏微分方程的工具箱函数求解

一.偏微分方程组的matlab求解语句
​	该命令用以求解以下的$PDE$方程式：

$c(x,t,u,\frac{\partial u }{\partial x})\frac{\partial u}{\partial t } = x^{-m}\frac{\partial (x^mf(x,t,u,\frac{\partial u}{\partial x}))}{\partial x }+s(x,t,u,\frac{\partial u }{\partial x})$

​	其中：$t \in [t_0,t_f],x \in [a,b]$。偏微分方程的初解：

$u(x,t_0) = v_0(x)$

​	边界条件为：

$p(x,t,u) +q(x,t)f(x,t,u,\frac{\partial u}{\partial x}) = 0$

​	下面介绍求解此类方程的函数用法：

$sol = pdepe(m,pdepe,icfun,bcfun,xmesh,tspan,options);$

​	$m:$对称参数。
​	$xmesh:$位置向量,$xmesh = [x_0,x_1,...x_N],x_0 = a,x_N = b$。
​	$tspan:$时间变量$t$的向量，$tspan = [t_0,t_1,...t_M],t_0 = t_0,t_M = t_f$。
​	$pdefun:$用户提供的$pde$函数文件。函数格式如下:

$[c,f,s] = pdefun(x,t,u,dudx);$

​					也就是说我们要自己设置相应的输出$c,f,s$，且它们都是行向量。
​	$icfun:$求解$u$的起始值，格式为$u = icfun(x)$。且$u$是行向量。
​	$bcfun:$提供边界条件函数，格式：

$[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t);$

​					$pl,ql:$左边界$p$和$q$的行向量。$pr,qr:$右边界$p$和$q$的行向量。
​	$options:$求解器相关解法参数，见$odeset$。
​	$sol:$多维向量输出，$sol(:,:,i)$为$u_i$的输出，而$u_i(j,k) = sol(j,k,i)$表示在$t = tspan(j),x = xmesh(k)$时候			的$u_i$的值。
​	要获得特定位置和时间的解用以下命令：

$[uout,duoutdx] = pdeval(m,xmesh,ui,xout);$

​	$xmesh:[x_0,x_1,...x_N]$
​	$ui:sol(j,:,i),$第$i$个输出$u_i$在时间$t_j$处的解。
​	$uout:$在指定$t_f$下对应指定位置$xout$的值。
​	$duoutdx:$相对应的$\frac{du}{dx}$。
二.具体的用法
​	1.求解以下偏微分方程(解析解为 $u(x,t) = e^{-t}sin(\pi x)$)：

$\pi^2\frac{\partial u}{\partial t} = \frac{\partial^2u}{\partial x^2}$

​		其中$x \in[0,1]$,满足以下条件:

$u(x,0) = sin(\pi x)\\
u(0,t) = 0 \\
\pi e^{-t} + \frac{\partial u(1,t)}{\partial x} = 0$
​	$solve:$
​	改写以上偏微分方程到标准形式：

$\pi^2 \frac{\partial u}{\partial t} = x^0\frac{\partial}{\partial x}(x^0\frac{\partial u}{\partial x}) +0$

​	具体的实现代码如下:
function first
    %计算从t:0~3的值
    x = linspace(0,1,20);
    t = linspace(0,3,60);
    subplot(121);
    sol = pdepe(0,@firstPdefun,@firstIcfun,@firstBcfun,x,t);
    u = surf(x,t,sol(:,:,1));
    title('微分方程数值解');
    xlabel('x');
    ylabel('t');
    zlabel('u')
    subplot(122);
    [X,T] = meshgrid(x,t);
    U = exp(-T).*sin(pi*X);
    surf(X,T,U);
    title('微分方程解析解');
end

%方程段
function [c,f,s] = firstPdefun(x,t,u,dudx)
    c = pi^2;
    f = dudx;
    s = 0;
end
%起始值条件段
function u = firstIcfun(x)
    u = sin(pi*x);
end
%边界条件段
function [pl,ql,pr,qr] = firstBcfun(xl,ul,xr,ur,t);
    pl = ul;
    ql = 0;
    pr = pi*exp(-t);
    qr = 1;
end

对比一下发现几乎和解析解一摸一样！


2.求解以下偏微分方程的数值解：

$\frac{\partial u_1}{\partial t} = 0.024\frac{\partial^2u_1}{\partial x^2} - F(u_1-u_2)\\
\frac{\partial u_2}{\partial t} = 0.170\frac{\partial^2u_2}{\partial x^2} + F(u_1-u_2)\\
F(u_1-u_2) = e^{5.73(u_1-u_2)} - e^{-11.46(u_1-u_2)}$

初值条件:

$u_1(x,0) = 1\\
u_2(x,0) = 0$

边值条件:

$\frac{\partial u_1(0,t)}{\partial x} = 0\\
u_2(0,t) = 0\\
u_1(1,t) = 1\\
\frac{\partial u_2(1,t)}{\partial x} = 0\\$

$Solve:$
​	化简上面的偏微分方程为标准形式：

$\begin{pmatrix} 1\\1 \end{pmatrix}_.*\frac{\partial}{\partial t}\begin{pmatrix} u_1\\u_2 \end{pmatrix} = \frac{\partial}{\partial x}\begin{pmatrix} 0.024\frac{\partial u_1}{\partial x}\\0.170\frac{\partial u_2}{\partial x} \end{pmatrix} +\begin{pmatrix} -F(u_1-u_2)\\F(u_1-u_2) \end{pmatrix}$

​	化简左边界条件也有：

$\begin{pmatrix} 0\\u_2 \end{pmatrix}+\begin{pmatrix} 1\\0 \end{pmatrix}_.* \begin{pmatrix} 0.024\frac{\partial u_1}{\partial x}\\0.170\frac{\partial u_2}{\partial x} \end{pmatrix}=  \begin{pmatrix} 0\\0 \end{pmatrix}$

​	化简右边界条件有：

$\begin{pmatrix} u_1-1\\0 \end{pmatrix}+\begin{pmatrix} 0\\1 \end{pmatrix}_.* \begin{pmatrix} 0.024\frac{\partial u_1}{\partial x}\\0.170\frac{\partial u_2}{\partial x} \end{pmatrix}=  \begin{pmatrix} 0\\0 \end{pmatrix}$

​	最后附上整个代码：
function second
xmesh  = linspace(0,1,20);
tspan = linspace(0,3,60);
sol = pdepe(0,@secondPdefun,@secondPdein,@secondPdebc,xmesh,tspan);
subplot(121);
surf(xmesh,tspan,sol(:,:,1));
xlabel('x');
ylabel('t');
zlabel('u(1)');
title('u(1)-x-t图像');
subplot(122);
surf(xmesh,tspan,sol(:,:,2));
xlabel('x');
ylabel('t');
zlabel('u(2)');
title('u(2)-x-t图像');
end

function [c,f,s] = secondPdefun(x,t,u,dudx)
c = [1 1]';
f = [0.024*dudx(1) 0.170*dudx(2)]';
y = u(1) - u(2);
stemp = exp(5.73*y) - exp(-11.46*y);
s = [-stemp stemp]';
end

function up = secondPdein(x,t,u,dudx)
up = [1 0]';
end

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


最后的结果是：

1 pdepe()函数的一般调用格式 [1]
sol=pdepe(m,@pdefun,@pdeinit,@pdebound,x,t)，其中pdefun是偏微分方的描述函数，标准形式为：
初始条件满足：
边界条件满足：
说明：
a,b表示上、下边界；
ua，ub是上、下边界的近似解；
pa,qa对应xa上计算的p,q值；
pb,qb对应xb上计算的p,q值；
2 一个实例
3 边界条件限定描述
4 结果展示
4求解代码
function sys=pdepeA()
%{
程序功能：
1、计算一般的偏微分方程组
2、pdepe()函数的一般调用格式是：
sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)，其中pdefun是偏微分方程的描述函数
有固定的格式。
3、pdeic是偏微分方程的初始条件，初始条件的描述为u(x,t0)=u0，
可以使用u0=pde(x);
4、pdebc是偏微分方程的边界条件，它的标准形式为：
p(x,t,u)+q(x,t,u).*f(x,t,u,ux)=0，则可以用[pa, qa, pb, qb]=pdebc(x,t,u,ux),a和b表示下边界和上边界。


参考链接：
https://jingyan.baidu.com/article/af9f5a2d15494e43140a45e3.html


%}
    x=0:0.05:1;
    t=0:0.05:2;
    m=0;
    sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);
    u1=sol(:,:,1);
    u2=sol(:,:,2);
    
    figure
    surf(x,t,u1)
    title('u1(x,t)')
    xlabel('Distance x')
    ylabel('Time t')
    
    figure
    surf(x,t,u2)
    title('u2(x,t)')
    xlabel('Distance x')
    ylabel('Time t')

    sys=[];

end
%方程形式
function [c, f, s]=pdefun(x, t, u, ux)
    c=[1;1] ;
    y=u(1)-u(2);
    F=exp(5.73*y)-exp(-11.46*y);
    s=F*[-1;1];
    f=[0.024*ux(1) ;0.017*ux(2) ];

end
%初始条件
function [u0]=pdeic(x)
    u0=[1 ;0];
    

end

%边界条件
function [pa, qa, pb, qb]=pdebc(xa ,ua ,xb ,ub, t)
    pa=[0; ua(2)];
    qa=[1; 0];
    pb=[ub(1)-1; 0];
    qb=[0; 1];

end
参考
^链接 https://jingyan.baidu.com/article/af9f5a2d15494e43140a45e3.html

题目写的看不懂。我也瞎答。
写了个MATLAB的小程序，用特征线法求解 偏微分方程组。
[XX,TT]=meshgrid(0:0.4:4,0:0.1:1);
N=size(XX,2);
T=size(XX,1);
u1=zeros(T,N);
u2=zeros(T,N);
dx=0.4;
dt=0.1;
ds=0.1;
u10=1:N;u20=(N:-1:1); % initial value
u10t=1:T;u20t=1:T; % value at a given node point
u1(1,:)=u10;
u2(1,:)=u20;
u1(:,1)=u10t;
u2(:,1)=u20t;
for tt=2:T % 特征线法求解线性偏微分方程组
for ii=2:N
u1(tt,ii)=(-400*u1(tt-1,ii-1)+400*u2(tt-1,ii-1))*ds;
u2(tt,ii)=(0.35*u1(tt-1,ii)-0.35*u2(tt-1,ii))*dt;
end
end
quiver(XX,TT,u1,u2)

	

       此教程说明如何使用 MATLAB 构造几种不同类型的微分方程并求解。MATLAB 提供了多种数值算法来求解各种微分方程：
初始值问题
边界值问题
时滞微分方程
偏微分方程
初始值问题
vanderpoldemo 是用于定义 van der Pol 方程的函数
type vanderpoldemo
function dydt = vanderpoldemo(t,y,Mu)
dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];
       该方程写作包含两个一阶常微分方程 (ODE) 的方程组。将针对参数 μ 的不同值计算这些方程。为了实现更快的积分，您应该根据 μ 的值选择合适的求解器。
       当 μ=1 时，任何 MATLAB ODE 求解器都能有效地求解 van der Pol 方程。ode45 求解器就是其中之一。该方程在域 [0,20] 中求解，初始条件为 y(0)=2 和 dydtt=0=0。
tspan = [0 20];
y0 = [2; 0];
Mu = 1;
ode = @(t,y) vanderpoldemo(t,y,Mu);
[t,y] = ode45(ode, tspan, y0);% Plot solution
plot(t,y(:,1))
xlabel('t')
ylabel('solution y')
title('van der Pol Equation, \mu = 1')
       对于较大的 μ，问题将变为刚性。此标签表示拒绝使用普通方法计算的问题。这种情况下，要实现快速积分，需要使用特殊的数值方法。ode15s、ode23s、ode23t 和 ode23tb 函数可有效地求解刚性问题。
      当 μ=1000 时，van der Pol 方程的求解使用 ode15s，初始条件相同。您需要将时间范围大幅度延长到 [0,3000] 才能看到解的周期性变化。
tspan = [0, 3000];
y0 = [2; 0];
Mu = 1000;
ode = @(t,y) vanderpoldemo(t,y,Mu);
[t,y] = ode15s(ode, tspan, y0);
plot(t,y(:,1))
title('van der Pol Equation, \mu = 1000')
axis([0 3000 -3 3])
xlabel('t')
ylabel('solution y')
边界值问题
bvp4c 和 bvp5c 可以求解常微分方程的边界值问题。
      示例函数 twoode 将一个微分方程写作包含两个一阶 ODE 的方程组。此微分方程为
type twoode
function dydx = twoode(x,y)
%TWOODE  Evaluate the differential equations for TWOBVP.
dydx = [ y(2); -abs(y(1)) ];
函数 twobc 求解该问题的边界条件为：y(0)=0 和 y(4)=−2。
type twobc
function res = twobc(ya,yb)
res = [ ya(1); yb(1) + 2 ];
       在调用 bvp4c 之前，您必须为要在网格中表示的解提供一个猜想值。然后，求解器就像对解进行平滑处理一样修改网格。
bvpinit 函数以您可以传递给求解器 bvp4c 的形式设定初始猜想值。对于 [0 1 2 3 4] 的网格以及 y(x)=1 和 y'(x)=0 的常量猜想值，对 bvpinit 的调用为：
solinit = bvpinit([0 1 2 3 4],[1; 0]);
       利用这个初始猜想值，您可以使用 bvp4c 对该问题求解。使用 deval 计算 bvp4c 在某些点返回的解，然后绘制结果值。
sol = bvp4c(@twoode, @twobc, solinit);
xint = linspace(0, 4, 50);
yint = deval(sol, xint);
plot(xint, yint(1,:));
xlabel('x')
ylabel('solution y')
hold on
      此特定的边界值问题实际上有两种解。通过将初始猜想值更改为 y(x)=−1 和 y'(x)=0，可以求出另一个解。
solinit = bvpinit([0 1 2 3 4],[-1; 0]);
sol = bvp4c(@twoode,@twobc,solinit);
xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:));
legend('Solution 1','Solution 2')
hold off
时滞微分方程
     dde23、ddesd 和 ddensd 可以求解具有各种时滞的时滞微分方程。示例 ddex1、ddex2、ddex3、ddex4 和 ddex5 构成了这些求解器的迷你使用教程。
ddex1 示例说明如何求解微分方程组
y′1(t)=y1(t−1)y′2(t)=y1(t−1)+y2(t−0.2)y′3(t)=y2(t).
您可以使用匿名函数表示这些方程
ddex1fun = @(t,y,Z) [Z(1,1); Z(1,1)+Z(2,2); y(2)];
问题的历史解(t≤0 时)固定不变：
y1(t)=1y2(t)=1y3(t)=1.
您可以将历史解表示为由 1 组成的向量。
ddex1hist = ones(3,1);
采用二元素向量表示方程组中的时滞。
lags = [1 0.2];
      将函数、时滞、历史解和积分区间 [0,5] 作为输入传递给求解器。求解器在整个积分区间生成适合绘图的连续解。
sol = dde23(ddex1fun, lags, ddex1hist, [0 5]);
plot(sol.x,sol.y);
title({'An example of Wille and Baker', 'DDE with Constant Delays'});
xlabel('time t');
ylabel('solution y');
legend('y_1','y_2','y_3','Location','NorthWest');
偏微分方程
     pdepe 使用一个空间变量和时间对偏微分方程求解。示例 pdex1、pdex2、pdex3、pdex4 和 pdex5 构成了 pdepe 的迷你使用教程。
此示例问题使用函数 pdex1pde、pdex1ic 和 pdex1bc。
pdex1pde 定义微分方程
type pdex1pde
function [c,f,s] = pdex1pde(x,t,u,DuDx)
%PDEX1PDE  Evaluate the differential equations components for the PDEX1 problem.
c = pi^2;
f = DuDx;
s = 0;
pdex1ic 设置初始条件
type pdex1ic
function u0 = pdex1ic(x)
%PDEX1IC  Evaluate the initial conditions for the problem coded in PDEX1.
u0 = sin(pi*x);
pdex1bc 设置边界条件
u(0,t)=0,
πe−t+∂∂xu(1,t)=0.
type pdex1bc
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
%PDEX1BC  Evaluate the boundary conditions for the problem coded in PDEX1.
pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;
pdepe 需要提供空间离散 x 和时间向量 t(您要获取解快照的时间点)。使用包含 20 个节点的网格求解此问题，并请求五个 t 值的解。提取解的第一个分量并绘图。
x = linspace(0,1,20);
t = [0 0.5 1 1.5 2];
sol = pdepe(0,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
u1 = sol(:,:,1);
surf(x,t,u1);
xlabel('x');
ylabel('t');
zlabel('u');
关注公众号：MATLAB基于模型的设计 (ID：xaxymaker) ，每天推送MATLAB学习最常见的问题，每天进步一点点，业精于勤荒于嬉。
可保存后扫码关注哦！