matlab求解偏微分方程的范例

以下是一个使用MATLAB求解偏微分方程的范例：

考虑二维热传导方程：

∂u/∂t = κ(∂^2u/∂x^2 + ∂^2u/∂y^2)

其中，u(x,y,t)表示温度分布，κ是热传导系数。

假设在一个矩形区域内，初始时刻温度分布为：

u(x,y,0) = sin(πx)sin(πy)

并且在边界上设定了绝热的边界条件，即：

u(0,y,t) = u(1,y,t) = u(x,0,t) = u(x,1,t) = 0

现在需要求解该方程在一定时间范围内的温度分布。

解决此问题的MATLAB代码如下：

% 定义参数
L = 1;              % 区域长度
T = 0.1;            % 时间范围
nx = 50;            % x方向网格点数
ny = 50;            % y方向网格点数
nt = 1000;          % 时间步数
dx = L/(nx-1);      % x方向网格步长
dy = L/(ny-1);      % y方向网格步长
dt = T/nt;          % 时间步长
x = linspace(0, L, nx);
y = linspace(0, L, ny);
[X, Y] = meshgrid(x, y);
u = zeros(nx, ny);

% 设置初始条件
u(:, :) = sin(pi*X).*sin(pi*Y);    % u(x,y,0) = sin(pi*x)*sin(pi*y)

% 设置边界条件
u(1, :) = 0;        % u(0,y,t) = 0
u(nx, :) = 0;       % u(1,y,t) = 0
u(:, 1) = 0;        % u(x,0,t) = 0
u(:, ny) = 0;       % u(x,1,t) = 0

% 迭代求解
for n = 1:nt
    un = u;
    % 采用有限差分法求解
    u(2:nx-1, 2:ny-1) = un(2:nx-1, 2:ny-1) + ...
        dt*(1/(dx^2)*(un(3:nx, 2:ny-1) - 2*un(2:nx-1, 2:ny-1) + un(1:nx-2, 2:ny-1)) + ...
        1/(dy^2)*(un(2:nx-1, 3:ny) - 2*un(2:nx-1, 2:ny-1) + un(2:nx-1, 1:ny-2)));
end

% 绘制温度分布图
surf(X, Y, u)
xlabel('x')
ylabel('y')
zlabel('Temperature')

运行这段代码，就可以得到在一定时间范围内的温度分布图。