I need help with Predictor-Corrector Scheme for PDE. I use 2d matrices instead of 3d matrices, but something is wrong and I don't understand what.

This is my problem.

∂u/∂t = 4∆u + e^t * cos(πx/2) * sin(πy/4) x ∈ (0, 1), y ∈ (0, 2), t ∈ (0, T],

∂u/∂x | x=0 = y,

u | x=1 = y + 1,

u | y=0 = 1,

∂u/∂y | y=2 = x,

u | t=0 = x*y + 1

% CODE

% Parameters

Lx = 1.0; Ly = 2.0;

Nx = 50; Ny = 100; % Number of steps in x and y directions

% Coefficient

D = 4;

% Time parameters

T = 0.0001;

Nt = 1000; % Number of time steps

dt = T / Nt;

% Grid

dx = Lx / (Nx - 0.5);

x = linspace(-dx / 2, Lx, Nx + 1);

dy = Ly / (Ny - 0.5);

y = linspace(Ly + dy / 2, 0, Ny + 1);

t = linspace(0, T, Nt + 1);

% Initial condition

[X, Y] = meshgrid(x, y);

u = X .* Y + 1;

% Inhomogeneity function

f = @(x, y, t) exp(t) .* cos(pi * x / 2) .* sin(pi * y / 4);

% Thomas algorithm for tridiagonal system

function x = thomas_algorithm(a, b, c, d)

n = length(d);

cp = zeros(n - 1, 1);

dp = zeros(n, 1);

% Forward sweep

cp(1) = c(1) / b(1);

dp(1) = d(1) / b(1);

for i = 2:n

m = 1.0 / (b(i) - a(i - 1) * cp(i - 1));

if i < n

cp(i) = c(i) * m;

end

dp(i) = (d(i) - a(i - 1) * dp(i - 1)) * m;

end

% Backward substitution

x = zeros(n, 1);

x(n) = dp(n);

for i = n - 1:-1:1

x(i) = dp(i) - cp(i) * x(i + 1);

end

end

% Numerical solution scheme

for n = 1:Nt

u_old = u;

% Predictor step in x direction

for j = 2:Ny

a = -D * ones(Nx - 1, 1) / (dx^2);

b = (1 / (dt / 2) + 2 * D / (dx^2)) * ones(Nx - 1, 1);

c = -D * ones(Nx - 1, 1) / (dx^2);

d = u(2:Nx, j) / (dt / 2) + f(x(2:Nx), y(j), t(n));

u(2:Nx, j) = thomas_algorithm(a, b, c, d);

end

% Predictor step in y direction

for i = 2:Nx

a = -D * ones(Ny - 1, 1) / (dy^2);

b = (1 / (dt / 2) + 2 * D / (dy^2)) * ones(Ny - 1, 1);

c = -D * ones(Ny - 1, 1) / (dy^2);

d = u(i, 2:Ny) / (dt / 2) + f(x(i), y(2:Ny), t(n));

u(i, 2:Ny) = thomas_algorithm(a, b, c, d);

end

u_old1 = u;

% Corrector step

for i = 2:Nx

for j = 2:Ny

u(i, j) = u_old(i, j) + dt * (...

D * (u_old1(i - 1, j) - 2 * u_old1(i, j) + u_old1(i + 1, j)) / (dx^2) + ...

D * (u_old1(i, j - 1) - 2 * u_old1(i, j) + u_old1(i, j + 1)) / (dy^2) + ...

f(x(i), y(j), t(n) + dt / 2) ...

);

end

end

% Apply boundary conditions

u(:, 1) = 1; % y = 0

u(:, end) = u(:, end - 1) + dy * x; % Neumann y = 2

u(1, :) = u(2, :) - dx * y; % Neumann x = 0

u(end, :) = y + 1; % x = 1

end

Answer 1

[Aastha]

Hi Jane,

I ran the code and encountered some errors. I have mentioned the errors and suggestions to resolve them below:

In the line:

d = u(2:Nx, j) / (dt / 2) + f(x(2:Nx), y(j), t(n));

There was an error that I diagnosed using a "try-catch" block. It appears that the variable "u" is of size 100x51, but the loop iterator “j” is going from 2 to “Ny” which is causing an out-of-bounds error since “Ny” is greater than 51. You can either reduce the value of “Ny” to 50 (since you are using “j+1” as an index for “u” in your code) or increase the size of “u” based on your application's needs.

Next, in this line:

u(:, end) = u(:, end - 1) + dy * x; % Neumann y = 2;

there is a size mismatch between the left-hand side (LHS) and the right-hand side (RHS). The term “u(:, end - 1)” is of size 51x1, while “dy * x” is of size 1x51. Adding these results in a 51x51 matrix, which doesn't match the LHS size of 51x1. You can resolve this by using the “transpose” function in MATLAB to make “dy * x” a 51x1 matrix with the code snippet mentioned below:

transpose(dy * x);

This should fix the errors in the code.

You may refer to MathWorks documentation of “transpose” and “try-catch” function for more information. Following are the links for the same:

transpose:  https://www.mathworks.com/help/matlab/ref/double.transpose.html

try-catch: https://www.mathworks.com/help/matlab/ref/try.html

I hope this helps!