MATLAB偏微分方程数值解: 求解一维偏微分方程

下面结合一个简单的实例介绍一维PDE的求解。

【例1】求解下面的PDE问题：

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

其中， $0\leq x\leq1$ ， $t\geq0$ 。

当t=0时，解满足初始条件

\[u(x,0)=\operatorname{sin}{}\pi x\]

当x=0和x=1时，解满足下面的边界条件

\[\begin{matrix} u(0,t)=0 \\ \pi{e}^{-t}+\frac{\partial u}{\partial x}(1,t)=0 \end{matrix}\]

按照下面的步骤求解此方程。

1．重写PDE

按照方程(1-1)的形式重写PDE，即

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

参数m=0，项

\[c(x,t,u,\frac{\partial u}{\partial x})={\pi}^{2}\]

\[f(x,t,u,\frac{\partial u}{\partial x})=\frac{\partial u}{\partial x}\]

\[s(x,t,u,\frac{\partial u}{\partial x})=0\]

下面把PDE问题写成程序，用函数pdexlpde表示。

code.matlab

function [c,f,s] = pdex1pde(x,t,u,DuDx)
c = pi^2;
f = DuDx;
s = 0;

编写初始条件的函数pdexlic，如下所示：

code.matlab

function u0 = pdex1ic(x)
u0 = sin(pi*x);

在MATLAB中编写边界条件的函数是pdexlbc，如下所示：

code.matlab

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;

其中，pl和ql对应于左侧的边界条件（x=0），pr和qr对应于右侧的边界条件（x=1）。

点(t,x)为希望pdepe函数求解的点。用向量t和x指定点。这两个向量在求解过程中起不同的作用。在计算时间和精度时与向量x的长度密切相关，而与向量t中的值则不是很相关。

本例中在[0,1]中取20个等间隔的点，在[0,2]中取5个等间隔的时间点，形成计算网格。下面计算网格节点上的解。

code.matlab

>> x = linspace(0,1,20);
>> t = linspace(0,2,5);

下面调用pdepe函数，m=0，使用函数pdex1pde、pdexlic、pdexlbc和x、t定义的计算网格。pdepe函数将数值解返回到三维数组sol中，其中，sol(i,j,k)近似表示解中的第k个组分uk，为t(i)和x(j)处的解。

code.matlab

>> m = 0;
>> sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);

本例使用函数句柄@来传递pdexlpde、pdexlic和pdexlbc函数。

下面提取和显示第一个解组分。在本例中，解u只有一个组分，但是为了达到演示的目的，本例从三维数组中将它“提取”出来。在命令窗口中输入下面的命令行

code.matlab

>> u = sol(:,:,1);
>>
>> surf(x,t,u)
>> title('方程的解')
>> xlabel('距离x')
>>ylabel('时间t')

结果如图1-1所示。

\[\]

图1-1 方程的解