issue in approximation of wave speed of propagation

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

[am]

@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

[Torsten]

The solution in the article is radially symmetric. So for each time tj, you only need to compute the radius rj for which the solution becomes > 0. The slope of this curve will give you the speed of propagation.

[am]

@Torsten thanks for your remark but i know that for Fisher kpp equation the speed of propagation is larger or equal to2 how can i show that numerically? In addition in one dimension the formula i give above in the first comment is true

[am]

% Parameters
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
c_range = linspace(0, 2, 5); % Narrower range of wave speeds to explore
nx = 51;            % Number of grid points in x-direction
ny = 51;            % Number of grid points in y-direction
nt = 5000;          % Number of time steps

% Grid spacing and time step
dx = Lx / (nx - 1);
dy = Ly / (ny - 1);
dt = C * dx;

% 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 array to store wave profiles
u_profiles = zeros(length(c_range), length(x));

% Loop over wave speeds
for i = 1:length(c_range)
    c = c_range(i);
    
    % Define time domain based on wave speed
    T = max(Lx, Ly) / c;
    timesteps = round(T / dt);
    
    % Initialize solution matrix
    u = zeros(1, length(x));

    % Loop over time steps
    for t = 1:timesteps
        % Update solution (if necessary)
    end

    % Store wave profile
    u_profiles(i, :) = u;
end

% Plot wave profiles
figure;
for j = 1:length(c_range)
    plot(x, u_profiles(j, :), 'DisplayName', sprintf('c = %.2f', c_range(j)));
    hold on;
end
xlabel('x');
ylabel('u(x)');
legend('show');

@Torsten i think that i need just to test values of wave sped to see if ithere is travelling wave solution sor not isuggest the code above :however im not sure if its a good way ti do that ,can you give me your opinion ,