Solving Partial Differential Equations in MATLAB

Solving Partial Differential Equations in MATLAB

Partial Differential Equations (PDEs) describe physical phenomena such as heat conduction, fluid dynamics, wave propagation, and electromagnetic fields. In MATLAB, you can solve PDEs through three primary approaches:

• The Built-in pdepe Solver: Solves initial-boundary value problems for systems of 1-D parabolic and elliptic PDEs.
• The Partial Differential Equation Toolbox: Solves 2-D and 3-D problems using the Finite Element Method (FEM) for structural mechanics, thermal analysis, electromagnetics, and general PDEs.
• Custom Numerical Methods: Solves custom discretized PDE domains using the Finite Difference Method (FDM) or Finite Volume Method (FVM).

1. Using pdepe (1-D PDEs)

The pdepe function solves systems of 1-D parabolic and elliptic PDEs across a spatial interval a ≤ x ≤ b and time interval t ≥ t₀.

Standard Mathematical Form

Your PDE must be written in the following standard form:

c(x, t, u, ∂u/∂x) · ∂u/∂t = x-m · ∂/∂x [ xm · f(x, t, u, ∂u/∂x) ] + s(x, t, u, ∂u/∂x)

Where:

• m: Coordinate symmetry parameter (m = 0 for Cartesian/slab, m = 1 for cylindrical, m = 2 for spherical).
• c: Capacity coefficient matrix or scalar.
• f: Flux term containing spatial derivative terms ∂u/∂x.
• s: Source term.
• u: Unknown state variable vector or scalar.

Standard Boundary Condition Form

Boundary conditions at the left boundary x = a and right boundary x = b must satisfy:

p(x, t, u) + q(x, t) · f(x, t, u, ∂u/∂x) = 0

Step-by-Step Workflow for pdepe

1. Define PDE coefficients in a function:
   MATLABCopy

   function [c, f, s] = mypde(x, t, u, dudx)
       c = 1;
       f = dudx; % Diffusion flux term
       s = 0;    % Source term
   end

2. Define initial conditions:
   MATLABCopy

   function u0 = myic(x)
       u0 = sin(pi * x); % Initial state at t = 0
   end

3. Define boundary conditions:
   MATLABCopy

   function [pl, ql, pr, qr] = mybc(xl, ul, xr, ur, t)
       % Left boundary at x = 0 (Dirichlet: u = 0)
       pl = ul;
       ql = 0;
       
       % Right boundary at x = 1 (Dirichlet: u = 0)
       pr = ur;
       qr = 0;
   end

4. Create spatial and time meshes:
   MATLABCopy

   x = linspace(0, 1, 50);   % 50 spatial points
   t = linspace(0, 0.5, 100); % 100 time evaluation points

5. Solve using pdepe:
   MATLABCopy

   m = 0; % Slab geometry
   sol = pdepe(m, @mypde, @myic, @mybc, x, t);

6. Extract the solution:
   MATLABCopy

   u = sol(:,:,1);
   % Solution matrix of size [length(t) x length(x)]

7. Plot the results:
   MATLABCopy

   figure;
   surf(x, t, u);
   xlabel('Position x');
   ylabel('Time t');
   zlabel('Solution u(x,t)');
   title('1-D Heat Diffusion Solution');
   shading interp;
   colorbar;

Complete 1-D Heat Equation Example (heat_equation_1d.m)

function solve_1d_heat()
    % Parameters
    m = 0; % 0 = slab/Cartesian
    x = linspace(0, 1, 40);
    t = linspace(0, 0.2, 50);

    % Solve PDE
    sol = pdepe(m, @pde_def, @pde_ic, @pde_bc, x, t);
    u = sol(:,:,1);

    % Visualization
    figure('Color', 'w');
    subplot(1,2,1);
    surf(x, t, u);
    xlabel('Distance x (m)');
    ylabel('Time t (s)');
    zlabel('Temperature u (^{circ}C)');
    title('Surface Plot of Temperature Profile');
    shading interp;
    colormap jet;
    colorbar;

    subplot(1,2,2);
    plot(x, u(1,:), 'k--', 'LineWidth', 1.5); hold on;
    plot(x, u(10,:), 'b-', 'LineWidth', 1.5);
    plot(x, u(25,:), 'g-', 'LineWidth', 1.5);
    plot(x, u(end,:), 'r-', 'LineWidth', 1.5);
    xlabel('Distance x (m)');
    ylabel('Temperature (^{circ}C)');
    title('Temperature at Different Time Steps');
    legend('t = 0 s', 't = 0.04 s', 't = 0.10 s', 't = 0.20 s', 'Location', 'NorthEast');
    grid on;
end

% PDE coefficient definitions: c * du/dt = d/dx(f) + s
function [c, f, s] = pde_def(x, t, u, dudx)
    alpha = 1.0; % Thermal diffusivity
    c = 1;
    f = alpha * dudx;
    s = 0;
end

% Initial condition u(x, 0)
function u0 = pde_ic(x)
    u0 = sin(pi * x) + 0.5 * sin(3 * pi * x);
end

% Boundary conditions: p + q * f = 0
function [pl, ql, pr, qr] = pde_bc(xl, ul, xr, ur, t)
    % Left side fixed at u = 0
    pl = ul - 0;
    ql = 0;
    
    % Right side fixed at u = 0
    pr = ur - 0;
    qr = 0;
end

2. Using Partial Differential Equation Toolbox (2-D and 3-D)

The PDE Toolbox uses the Finite Element Method (FEM) to solve scalar PDEs and systems of PDEs on 2-D and 3-D geometric models.

Standard PDE Models in the Toolbox

• Elliptic Equation (Steady-state heat, electrostatics): -∇ · (c ∇u) + au = f
• Parabolic Equation (Transient heat conduction): d · (∂u/∂t) - ∇ · (c ∇u) + au = f
• Hyperbolic Equation (Structural vibrations, wave propagation): m · (∂2u/∂t2) - ∇ · (c ∇u) + au = f

Step-by-Step Workflow for PDE Toolbox

1. Create the PDE model container:
   MATLABCopy

   model = createpde(); % Generic PDE model
   % Or physics-specific: model = createpde('thermal', 'transient-steady');

2. Define or import geometry:
   MATLABCopy

   % Built-in 2D shape
   geometryFromEdges(model, @lshapeg);
   
   % Or import 3D CAD geometries (STL or STEP format)
   % importGeometry(model, 'bracket.stl');

3. Specify PDE coefficients:
   MATLABCopy

   specifyCoefficients(model, 'm', 0, 'd', 0, 'c', 1.5, 'a', 0, 'f', 20);

4. Apply boundary conditions:
   MATLABCopy

   % Dirichlet boundary condition: fixed temperature u = 0
   applyBoundaryCondition(model, 'dirichlet', 'Edge', [1, 2, 3], 'u', 0);
   
   % Neumann boundary condition: insulated wall (flux g = 0)
   applyBoundaryCondition(model, 'neumann', 'Edge', [4, 5, 6], 'g', 0, 'q', 0);

5. Generate the finite element mesh:
   MATLABCopy

   mesh = generateMesh(model, 'Hmax', 0.08, 'GeometricOrder', 'quadratic');

6. Solve the PDE model:
   MATLABCopy

   results = solvepde(model); % Stationary solver
   % For time-dependent problems: results = solvepde(model, tlist);

7. Visualize the solution:
   MATLABCopy

   figure;
   pdeplot(model, 'XYData', results.NodalSolution, 'Contour', 'on');
   title('2-D Temperature Distribution');
   colormap jet;
   colorbar;

Complete 2-D Steady-State Heat Equation Example

% Step 1: Create PDE Model
model = createpde();

% Step 2: Define Geometry (L-shaped domain)
geometryFromEdges(model, @lshapeg);

% Step 3: Specify PDE Coefficients (-div(c*grad(u)) + a*u = f)
thermal_conductivity = 1.5; % c
heat_generation = 20.0;     % f
specifyCoefficients(model, 'm', 0, 'd', 0, ...
                           'c', thermal_conductivity, ...
                           'a', 0, ...
                           'f', heat_generation);

% Step 4: Apply Boundary Conditions
% Edges 1, 2, 3 held at 0 degrees C (Dirichlet)
applyBoundaryCondition(model, 'dirichlet', 'Edge', [1, 2, 3], 'u', 0);

% Remaining edges 4, 5, 6 are insulated (Zero Neumann flux)
applyBoundaryCondition(model, 'neumann', 'Edge', [4, 5, 6], 'g', 0, 'q', 0);

% Step 5: Generate Finite Element Mesh
mesh = generateMesh(model, 'Hmax', 0.08);

% Step 6: Solve the Model
results = solvepde(model);
u_solution = results.NodalSolution;

% Step 7: Plot Results
figure('Color', 'w');
subplot(1,2,1);
pdeplot(model, 'XYData', u_solution, 'Contour', 'on');
title('Temperature Distribution (^{circ}C)');
xlabel('X Coordinate (m)');
ylabel('Y Coordinate (m)');
colormap jet;
colorbar;

subplot(1,2,2);
pdemesh(model);
title(['Triangular FEM Mesh (' num2str(size(mesh.Nodes, 2)) ' Nodes)']);
xlabel('X Coordinate (m)');
ylabel('Y Coordinate (m)');
axis equal;

3. Solving PDEs with Finite Difference Method (Direct Scripting)

For custom boundary conditions or structured Cartesian grids, you can discretize PDEs directly using the Finite Difference Method (FDM) and solve the linear system A * u = b using MATLAB sparse matrices.

Example: 2-D Poisson Equation (∇²u = f(x,y))

% Grid Setup
N = 60; % Grid points per axis
L = 1.0;
h = L / (N - 1);
[X, Y] = meshgrid(linspace(0, L, N), linspace(0, L, N));

% Source term: Gaussian distribution at center
f = 100 * exp(-((X - 0.5).^2 + (Y - 0.5).^2) / 0.02);

% Assemble 1D second derivative operator
main_diag = -2 * ones(N, 1);
off_diag  = ones(N, 1);
D1 = spdiags([off_diag, main_diag, off_diag], [-1, 0, 1], N, N);

% Assemble 2D 5-point Laplacian using Kronecker tensor product
I = speye(N);
L2D = (kron(I, D1) + kron(D1, I)) / (h^2);

% Right hand side vector
rhs = reshape(f, N*N, 1);

% Apply Dirichlet boundary conditions (u = 0 on borders)
is_boundary = false(N, N);
is_boundary(1, :) = true;
is_boundary(N, :) = true;
is_boundary(:, 1) = true;
is_boundary(:, N) = true;
b_indices = find(is_boundary(:));

% Modify matrix rows for boundary conditions
for idx = b_indices'
    L2D(idx, :) = 0;
    L2D(idx, idx) = 1;
    rhs(idx) = 0;
end

% Solve linear system with sparse backslash
u_vec = L2D  rhs;
U = reshape(u_vec, N, N);

% Visualization
figure('Color', 'w');
surf(X, Y, U);
shading interp;
colormap turbo;
colorbar;
xlabel('X');
ylabel('Y');
zlabel('Potential u(x,y)');
title('2-D Poisson Equation Solved via Finite Difference Method');

4. Summary of Boundary Condition Formulations in MATLAB

5. Practical Tips for Solving PDEs in MATLAB

• Handling Discontinuities: When initial conditions have sharp gradients or step changes, refine your spatial grid near the step to prevent numerical oscillations.
• Mesh Convergence Checks: In 2D and 3D FEM models, verify solution accuracy by reducing the Hmax mesh parameter by half and confirming that peak temperature or stress values stabilize.
• Scalar Input in pdepe: The pdepe solver passes scalar values for x, t, u, and dudx. Ensure your coefficient subfunctions perform scalar arithmetic.
• Coordinate Symmetry Parameter: Set m = 0 for Cartesian slabs, m = 1 for cylindrical geometry (radial heat conduction in pipes), and m = 2 for spherical symmetry.