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
endThis 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 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 pde, predictor-corrector, 2d MATLAB Answers — New Questions