Explicit method for Allen-Cahn equation
12 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
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.
Réponse acceptée
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
Plus de réponses (0)
Voir également
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!