@Torsten thanks for your answer in my previous question ,your comments lmake me notice that its not necessarut the value 0.5 but its sufficient to have ;u(x_i,y,t_i)=u(x_j,y,t_j) , i also forget to tell you that im comparing my result with the resylts in the follox=wing article :Numerical Solution of the 1 + 2 Dimensional Fisher’s Equation by Finite Elements and the Galerkin Method by J. ROESSLER AND H. H~~SSNER
issue in approximation of wave speed of propagation
16 views (last 30 days)
Show older comments
0
im exploring Fisher KPP equation and diffrent numerical methods to study its behaviour specifically travelling waves and speed of propagation
here is the code in Matlab where iuse finite difference scheme:
C = 0.05; % Courant number
D = 1; % Diffusion coefficient
r = 1; % Reaction rate
Lx = 10; % Width of the domain
Ly = 10; % Height of the domain
nx = 51; % Number of grid points in x-direction
ny = 51; % Number of grid points in y-direction
nt = 5000; % Number of time steps
dx = Lx / (nx-1); % Spacing of grid in x-direction
dy = Ly / (ny-1); % Spacing of grid in y-direction
dt = C * dx; % Time step size
% Initialize concentration profile
un = zeros(ny, nx);
x = linspace(-Lx/2, Lx/2, nx);
y = linspace(-Ly/2, Ly/2, ny);
[X, Y] = meshgrid(x, y);
% Set initial condition
initial_condition = (X.^2 + 2*Y.^2 < 0.25);
un = 0.5 * initial_condition;
% Initialize matrix to save solutions over time
U = zeros(nt, ny, nx);
U(1,:,:) = un;
% Loop over time steps to solve the PDE
for n = 1:nt
uc = un;
% Update concentration profile using finite difference method
for i = 2:nx-1
for j = 2:ny-1
un(j, i) = uc(j, i) + ...
dt * (uc(j-1, i) - 2 * uc(j, i) + uc(j+1, i))/dy^2 + ...
dt * (uc(j, i-1) - 2 * uc(j, i) + uc(j, i+1))/dx^2 + ...
dt * uc(j, i) * (1 - uc(j,i));
end
end
% Apply Dirichlet boundary conditions
un(1,:) = 0;
un(end,:) = 0;
un(:,1) = 0;
un(:,end) = 0;
% Store solution
U(n+1,:,:) = un;
end
% Initialize arrays to store pairs of (x, t)
x_values = [];
t_values = [];
% Loop over time steps
for ti = 1:nt
for tj = ti+1:nt % to avoid comparing the same time step twice
% Extract concentration profiles at times ti and tj
ui = squeeze(U(ti,:,:));
uj = squeeze(U(tj,:,:));
% Find pairs (xi, xj) where ui(xi, y) = uj(xj, y) for all y
[xi, xj] = find(all(abs(ui - uj) < 1e-2, 1)); % tolerance for numerical solution
% Store pairs (xi, ti) and (xj, tj)
for k = 1:length(xi)
x_values = [x_values; x(xi(k)), x(xj(k))];
t_values = [t_values; ti, tj];
end
end
end
% Plot pairs of (ti, tj) against (xi, xj)
figure;
plot(t_values, x_values, 'o');
xlabel('Temps');
ylabel('Position x');
title('Paires (ti, tj) où u(xi, y, ti) = u(xj, y, tj) pour tout y');
grid on;
i'm expecting that the slope graph equal or large than 2 but with this code i dont get that ? what is the issue? the slope graph her is the approximate speed whixh u calculte by :c=(x_j-x_i)/(t_j-t_i) such that u(x_i,y,t_i)=u(x_j,y,t_j)can you please tell me the issue and if you have another suggestion please tell me
4 Comments
Answers (0)
See Also
Categories
Find more on Fit Postprocessing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!