Explicit method for Allen-Cahn equation

12 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)

Categorie

Scopri di più su Get Started with Curve Fitting Toolbox in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by