Explicit method for Allen-Cahn equation

7 visualizzazioni (ultimi 30 giorni)
Erm
Erm il 1 Ott 2024
Commentato: Erm il 1 Ott 2024
The plot of the equation must start at x=-1 and end at x=1. but mu result did not show that?
clear all;
clc;
maxk = 1000;
T = 0.10;
n = 50;
L = 2; % Length of the spatial domain [−1, 1]
Nx = 400; % Number of spatial grid points
dx = L / (Nx - 1); % Spatial step size
dt = T/maxk;
T = 1; % Final time
Nt = round(T / dt); % Number of time steps
a = 0.0001;
r = a * dt / (dx * dx); % Diffusion factor for explicit scheme
% Initial condition
x = linspace(-1, 1, n+1);
u = zeros(n+1, maxk+1);
u(:,1) = x.^2 .* cos(pi * x);
% Implementation of the explicit method for Allen-Cahn equation
for t = 1:maxk
% Internal points
for i = 2:n
u(i, t+1) = u(i, t) + r * (u(i-1, t) - 2 * u(i, t) + u(i+1, t)) ...
+ dt * (5 * u(i, t)^3 - 5 * u(i, t));
end
% Periodic boundary conditions
u(1, t+1) = u(end-1, t+1); % Periodic condition for first point
u(end, t+1) = u(2, t+1); % Periodic condition for last point
end
% Plot results
figure; % Create a new figure
xx = linspace(-1, 1, 100);
t_values = [0, 0.2, 0.4, 0.6, 0.8]; % Time values to plot
plot(x, u(:,1), '-', x, u(:,round(maxk*0.2)), '-', x, u(:,round(maxk*0.4)), '-', x, u(:,round(maxk*0.6)), '-', x, u(:,end), '-');
xlabel('x');
ylabel('u(x,t)');
grid on;
legend('t = 0', 't = 0.2', 't = 0.4', 't = 0.6', 't = 0.8');
hold off;
  1 Commento
Torsten
Torsten il 1 Ott 2024
Did you read somewhere that you need a boundary condition on the gradients of u ? In my opinion, u(-1,t) = u(1,t) suffices to fix a solution.

Accedi per commentare.

Risposta accettata

Torsten
Torsten il 1 Ott 2024
xstart = -1.0;
xend = 1.0;
nx = 401;
x = linspace(xstart,xend,nx).';
dx = x(2)-x(1);
tstart = 0.0;
tend = 1.0;
nt = 10;
tspan = linspace(tstart,tend,nt);
u0 = x.^2.*cos(pi*x);
D = 1e-4;
[T,U] = ode15s(@(t,u)fun(t,u,D,nx,dx),tspan,u0);
plot(x,[U(1,:);U(end,:)])
function dudt = fun(t,u,D,nx,dx)
ufull = [u(end-1);u;u(2)];
dudt = D*(ufull(3:end)-2*ufull(2:end-1)+ufull(1:end-2))/dx^2-5*ufull(2:end-1).^3+5*ufull(2:end-1);
end
  4 Commenti
Torsten
Torsten il 1 Ott 2024
Modificato: Torsten il 1 Ott 2024
On your code why you wrote D = 1e-4
I prefer the scientific notation for small numbers. Think about how many 0's you had to write if D was 1e-10 instead of 1e-4 :-)
Erm
Erm il 1 Ott 2024
got it. thank you.

Accedi per commentare.

Più risposte (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by