Azzera filtri
Azzera filtri

traveling waves of fisher kpp equation

39 visualizzazioni (ultimi 30 giorni)
am
am il 17 Mag 2024 alle 7:34
Modificato: am il 17 Mag 2024 alle 14:36
hi,i how to calculate wave speed of solution of Fisher Kpp equation:$u_{t}+u_{xx}+u_{yy}=u(1-u)$ i write a coode in matlab in which first i do finite difference scheme and calculate speed
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 = 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^2 / D; % Time step size, adjusted for stability
% 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;
% Variables to track wave front
front_positions = [];
front_times = [];
% 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 * D * (uc(j-1, i) - 4 * uc(j, i) + uc(j+1, i) + ...
uc(j, i-1) + uc(j, i+1)) / (dx * dx) + ...
dt * r * 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;
% Determine the front positions at half maximum concentration
half_max = 0.5;
half_max_positions = find(abs(un - half_max) < 0.01);
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
% Calculate the speed of the wave
if length(front_positions) > 1
front_velocity = diff(front_positions) ./ diff(front_times);
% Calculate the average speed of the wave
average_speed = mean(front_velocity);
else
average_speed = 0; % Not enough data to compute speed
end
% Display the average speed of the wave
disp('Average speed of the wave:');
Average speed of the wave:
disp(average_speed);
NaN
i dont know if the problem in the way i calculate the speed please can you give me a hint?

Risposte (1)

Torsten
Torsten il 17 Mag 2024 alle 13:27
Modificato: Torsten il 17 Mag 2024 alle 13:47
Your equation reads
u_{t} = - (u_{xx}+u_{yy}) + u(1-u)
You miss the minus sign in your code (at least if the equation you wrote is correct).
Further you should save the history of your time stepping results. The way you coded the problem, you can only plot the results of the last two time steps from uc and un.
If you look at the results over time, I don't see a propagating wave. Only the shape of the initial bell changes over time.
  1 Commento
am
am il 17 Mag 2024 alle 14:36
Modificato: am il 17 Mag 2024 alle 14:36
thanks for your answer @Torsten the equation in the code is right ,the right equation is u_{t} = (u_{xx}+u_{yy}) + u(1-u) this is a Fisher Kpp and it has travelling waves for wave speed c>=2 where c indicates how fast the wave profile is moving through the space ,for one dimension c is the rate at which the position x of a fixed point on the wave profile (e.g., where 𝑢=0.5) changes with time.i think that for two dimension i would hacve two speeds one in direction x and the other in direction y,do you see how tocalculte it?
thanks

Accedi per commentare.

Prodotti


Release

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by