issue in calculation of wave speed

hi, i'm aiming to calculate wave speed of solution of "u_{t}=u_{xx}+u_{yy}+u-u^2" so i use first finite difference scheme but my issue is that wen i run the code i get error "unrecognized function'average speed' how to solve that? THANKS

D = 1;               % Diffusion coefficient
r = 1;               % Reaction rate
Lx = 100;             % Width of the domain
Ly = 100;             % Height of the domain

nx = 510;             % Number of grid points in x-direction
ny = 510;        % Parameters
     % Number of grid points in y-directionn
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

C = 0.05;            % Courant number
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);

% Initial concentration
un(:,:,:) = exp(-X.^2 - Y.^2 );
t = 0;

% Loop over time steps
for n = 1:nt
    uc = un;
    t = t + dt;

    for i = 2:nx-1
        for j = 2:ny-1
             un(j, i) = uc(j, i) +...
             dt *(uc(j-1, i) - 4 * uc(j, i) + uc(j+1, i) + ...
             uc(j, i-1)  + uc(j, i+1)) / (dy * dy) + ...
             dt * uc(j, i) * (1 - uc(j,i)) / (dy * dy);
            end
        
    end

    % Apply Dirichlet boundary conditions
    un(1,:,:) = 0;
    un(end,:,:) = 0;
    un(:,1,:) = 0;
    un(:,end,:) = 0;
    un(:,:,1) = 0;
    un(:,:,end) = 0;
    % Store front positions and times
    if ~isempty(half_max_positions)
        [front_rows, front_cols] = ind2sub(size(un), half_max_positions);
        front_positions = [front_positions; (front_cols - 1) * dx];
        front_times = [front_times; ones(size(front_cols)) * t];
    end
end

Unrecognized function or variable 'half_max_positions'.

% Calculate the speed of the wave
front_velocity = diff(front_positions) ./ diff(front_times);

% Calculate the average speed of the wave
average_speed = mean(front_velocity);

% Display the average speed of the wave
disp('Average speed of the wave:');
disp(average_speed);

[Walter Roberson]

You have a problem with half_max_positions

Answer 1

[Torsten]

The below code will at least give you the correct update of the solution.

You don't initialize "half_max_positions" ; that's why you get the error.

This is not how a stable time stepsize for explicit Euler is computed:

C = 0.05;            % Courant number
dt = C * dx          % Time step size

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-directionn
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

C = 0.05;            % Courant number
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);

% Initial concentration
un = exp(-X.^2 - Y.^2 );

% Initialize matrix to save solutions over time
U = zeros(nt,ny,nx);
U(1,:,:) = un;

t = 0;

% Loop over time steps
for n = 1:nt
    uc = un;
    t = t + dt;

    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(nt+1,:,:) = un;

    % Store front positions and times
    %if ~isempty(half_max_positions)
    %    [front_rows, front_cols] = ind2sub(size(un), half_max_positions);
    %    front_positions = [front_positions; (front_cols - 1) * dx];
    %    front_times = [front_times; ones(size(front_cols)) * t];
    %end
end
surf(X,Y,squeeze(U(nt+1,:,:)))

% Calculate the speed of the wave
%front_velocity = diff(front_positions) ./ diff(front_times);

% Calculate the average speed of the wave
%average_speed = mean(front_velocity);

% Display the average speed of the wave
%disp('Average speed of the wave:');
%disp(average_speed);

[am]

@Torsten Thank you how can i find values of positions x such that u(x,y,t)=0.5 and then dress a function x(t) (where x is the position where u(x,y,t)=0.5)?

[Torsten]

If you use

contourf(X,Y,squeeze(U(nt+1,:,:)),[0.5 0.5])

instead of

surf(X,Y,squeeze(U(nt+1,:,:)))

in the above code, you will see the isoline where u = 0.5.

I'm not sure what you want to extract here for x for the different time steps.

Maybe something like this:

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-directionn
nt = 500;            % Number of time steps
dx = Lx / (nx-1);    % Spacing of grid in x-direction
dy = Ly / (ny-1);    % Spacing of grid in y-direction

C = 0.05;            % Courant number
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);

% Initial concentration
un = exp(-X.^2 - Y.^2 );

% Initialize matrix to save solutions over time
U = zeros(nt,ny,nx);
U(1,:,:) = un;

t = 0;

% Loop over time steps
for n = 1:nt
    uc = un;
    t = t + dt;

    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
level = 0.3:0.1:0.7;
hold on
for i = 1:numel(level)
  xmax = nan(nt+1,1);
  for n = 1:nt+1
    M = contourf(X,Y,squeeze(U(n,:,:)),[level(i) level(i)]);
    if ~isempty(M)
      xmax(n) = max(abs(M(1,2:end)));
    end
  end
  plot(1:nt+1,xmax)
end
hold off

[am]

@Torstenthanks , the graph of x(t) where x and t are such that u(x,y,t)=0.5 it slope is an approximate value of speed of travelling waves of this equation ;my aim is to have the graph slope equal or large than 2 so i try to modify the initial condiition to :% Set initial condition

un(X == 0.5 & Y == 0.5) = 0.05;

un(X ~= 0.5 | Y ~= 0.5) = 0 i took this value from an article to be sure that it works ;but the result is worse what is not logical because here the graph solpe is an approximate spped of travelling waves of Fisherkpp which is know to be at least 2 so theere is something wrong in my code