ADI method for 2D reaction diffusion equation
hi , i want to use ADI mathod to fisherkpp equation and calcultae the integral between t=0 and t=10 however im not sure that i apply the method correctly ? so i want your feedbacks and your suggestion sto improve the code
here the code
% Parameters
Lx = 60; % Length of domain in x-direction
Ly = 1; % Length of domain in y-direction
Nx = 600; % Number of grid points in x-direction
Ny = 10; % Number of grid points in y-direction
final_time = 10; % Final time
D = 1; % Diffusion coefficient
r = 1; % Growth rate
% Discretization
dx = Lx / (Nx – 1);
dy = Ly / (Ny – 1);
dt = 0.006; % Time step
Nt = 5500; % Number of time steps
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
% Initial condition
[X, Y] = meshgrid(x, y);
u0 = (1/2 – 1/2 * tanh(1/(2*sqrt(6)) * X)).^2;
% ADI method
alpha = D * dt / dx^2;
beta = D * dt / dy^2;
gamma = r * dt;
u = u0;
u_new = u;
integral_u = 0;
for t = 1:Nt
% Implicit step in x-direction
for j = 2:Ny-1
% Construct tridiagonal matrix for x-direction
A = -alpha * ones(1, Nx-3);
B = (1 + 2 * alpha) * ones(1, Nx-2);
C = -alpha * ones(1, Nx-3);
D_vec = u(j, 2:end-1) + gamma * u(j, 2:end-1) .* (1 – u(j, 2:end-1));
% Solve tridiagonal system
u_new(j, 2:end-1) = tridiag_solve(A, B, C, D_vec);
end
% Implicit step in y-direction
for i = 2:Nx-1
% Construct tridiagonal matrix for y-direction
A = -beta * ones(1, Ny-3);
B = (1 + 2 * beta) * ones(1, Ny-2);
C = -beta * ones(1, Ny-3);
D_vec = u_new(2:end-1, i) + gamma * u_new(2:end-1, i) .* (1 – u_new(2:end-1, i));
% Solve tridiagonal system
u_new(2:end-1, i) = tridiag_solve(A, B, C, D_vec);
end
% Update integral of u over domain
if t > 1
integral_u = integral_u + trapz(y, trapz(x, u_new, 2), 1) * dt; % Trapezoidal numerical integration
end
u = u_new;
end
disp([‘Integral of u(x, y) over domain between t=0 and final time: ‘, num2str(integral_u)]);
function x = tridiag_solve(A, B, C, D)
% Solves a tridiagonal system of equations Ax = D
N = length(B);
c = zeros(N-1, 1);
d = zeros(N, 1);
x = zeros(N, 1);
c(1) = C(1) / B(1);
d(1) = D(1) / B(1);
for i = 2:N-1
denom = B(i) – A(i-1) * c(i-1);
c(i) = C(i) / denom;
d(i) = (D(i) – A(i-1) * d(i-1)) / denom;
end
d(N) = (D(N) – A(N-1) * d(N-1)) / (B(N) – A(N-1) * c(N-1));
x(N) = d(N);
for i = N-1:-1:1
x(i) = d(i) – c(i) * x(i+1);
end
endhi , i want to use ADI mathod to fisherkpp equation and calcultae the integral between t=0 and t=10 however im not sure that i apply the method correctly ? so i want your feedbacks and your suggestion sto improve the code
here the code
% Parameters
Lx = 60; % Length of domain in x-direction
Ly = 1; % Length of domain in y-direction
Nx = 600; % Number of grid points in x-direction
Ny = 10; % Number of grid points in y-direction
final_time = 10; % Final time
D = 1; % Diffusion coefficient
r = 1; % Growth rate
% Discretization
dx = Lx / (Nx – 1);
dy = Ly / (Ny – 1);
dt = 0.006; % Time step
Nt = 5500; % Number of time steps
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
% Initial condition
[X, Y] = meshgrid(x, y);
u0 = (1/2 – 1/2 * tanh(1/(2*sqrt(6)) * X)).^2;
% ADI method
alpha = D * dt / dx^2;
beta = D * dt / dy^2;
gamma = r * dt;
u = u0;
u_new = u;
integral_u = 0;
for t = 1:Nt
% Implicit step in x-direction
for j = 2:Ny-1
% Construct tridiagonal matrix for x-direction
A = -alpha * ones(1, Nx-3);
B = (1 + 2 * alpha) * ones(1, Nx-2);
C = -alpha * ones(1, Nx-3);
D_vec = u(j, 2:end-1) + gamma * u(j, 2:end-1) .* (1 – u(j, 2:end-1));
% Solve tridiagonal system
u_new(j, 2:end-1) = tridiag_solve(A, B, C, D_vec);
end
% Implicit step in y-direction
for i = 2:Nx-1
% Construct tridiagonal matrix for y-direction
A = -beta * ones(1, Ny-3);
B = (1 + 2 * beta) * ones(1, Ny-2);
C = -beta * ones(1, Ny-3);
D_vec = u_new(2:end-1, i) + gamma * u_new(2:end-1, i) .* (1 – u_new(2:end-1, i));
% Solve tridiagonal system
u_new(2:end-1, i) = tridiag_solve(A, B, C, D_vec);
end
% Update integral of u over domain
if t > 1
integral_u = integral_u + trapz(y, trapz(x, u_new, 2), 1) * dt; % Trapezoidal numerical integration
end
u = u_new;
end
disp([‘Integral of u(x, y) over domain between t=0 and final time: ‘, num2str(integral_u)]);
function x = tridiag_solve(A, B, C, D)
% Solves a tridiagonal system of equations Ax = D
N = length(B);
c = zeros(N-1, 1);
d = zeros(N, 1);
x = zeros(N, 1);
c(1) = C(1) / B(1);
d(1) = D(1) / B(1);
for i = 2:N-1
denom = B(i) – A(i-1) * c(i-1);
c(i) = C(i) / denom;
d(i) = (D(i) – A(i-1) * d(i-1)) / denom;
end
d(N) = (D(N) – A(N-1) * d(N-1)) / (B(N) – A(N-1) * c(N-1));
x(N) = d(N);
for i = N-1:-1:1
x(i) = d(i) – c(i) * x(i+1);
end
end hi , i want to use ADI mathod to fisherkpp equation and calcultae the integral between t=0 and t=10 however im not sure that i apply the method correctly ? so i want your feedbacks and your suggestion sto improve the code
here the code
% Parameters
Lx = 60; % Length of domain in x-direction
Ly = 1; % Length of domain in y-direction
Nx = 600; % Number of grid points in x-direction
Ny = 10; % Number of grid points in y-direction
final_time = 10; % Final time
D = 1; % Diffusion coefficient
r = 1; % Growth rate
% Discretization
dx = Lx / (Nx – 1);
dy = Ly / (Ny – 1);
dt = 0.006; % Time step
Nt = 5500; % Number of time steps
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
% Initial condition
[X, Y] = meshgrid(x, y);
u0 = (1/2 – 1/2 * tanh(1/(2*sqrt(6)) * X)).^2;
% ADI method
alpha = D * dt / dx^2;
beta = D * dt / dy^2;
gamma = r * dt;
u = u0;
u_new = u;
integral_u = 0;
for t = 1:Nt
% Implicit step in x-direction
for j = 2:Ny-1
% Construct tridiagonal matrix for x-direction
A = -alpha * ones(1, Nx-3);
B = (1 + 2 * alpha) * ones(1, Nx-2);
C = -alpha * ones(1, Nx-3);
D_vec = u(j, 2:end-1) + gamma * u(j, 2:end-1) .* (1 – u(j, 2:end-1));
% Solve tridiagonal system
u_new(j, 2:end-1) = tridiag_solve(A, B, C, D_vec);
end
% Implicit step in y-direction
for i = 2:Nx-1
% Construct tridiagonal matrix for y-direction
A = -beta * ones(1, Ny-3);
B = (1 + 2 * beta) * ones(1, Ny-2);
C = -beta * ones(1, Ny-3);
D_vec = u_new(2:end-1, i) + gamma * u_new(2:end-1, i) .* (1 – u_new(2:end-1, i));
% Solve tridiagonal system
u_new(2:end-1, i) = tridiag_solve(A, B, C, D_vec);
end
% Update integral of u over domain
if t > 1
integral_u = integral_u + trapz(y, trapz(x, u_new, 2), 1) * dt; % Trapezoidal numerical integration
end
u = u_new;
end
disp([‘Integral of u(x, y) over domain between t=0 and final time: ‘, num2str(integral_u)]);
function x = tridiag_solve(A, B, C, D)
% Solves a tridiagonal system of equations Ax = D
N = length(B);
c = zeros(N-1, 1);
d = zeros(N, 1);
x = zeros(N, 1);
c(1) = C(1) / B(1);
d(1) = D(1) / B(1);
for i = 2:N-1
denom = B(i) – A(i-1) * c(i-1);
c(i) = C(i) / denom;
d(i) = (D(i) – A(i-1) * d(i-1)) / denom;
end
d(N) = (D(N) – A(N-1) * d(N-1)) / (B(N) – A(N-1) * c(N-1));
x(N) = d(N);
for i = N-1:-1:1
x(i) = d(i) – c(i) * x(i+1);
end
end numerical integration, numerical analysis MATLAB Answers — New Questions