Explicit method for Allen-Cahn equation

12 vues (au cours des 30 derniers jours)
Erm
Erm le 1 Oct 2024
Commenté : Erm le 1 Oct 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 commentaire
Torsten
Torsten le 1 Oct 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.

Connectez-vous pour commenter.

Réponse acceptée

Torsten
Torsten le 1 Oct 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 commentaires
Torsten
Torsten le 1 Oct 2024
Modifié(e) : Torsten le 1 Oct 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 le 1 Oct 2024
got it. thank you.

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Get Started with Curve Fitting Toolbox dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by