traveling waves of fisher kpp equation
25 次查看(过去 30 天)
显示 更早的评论
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:');
disp(average_speed);
i dont know if the problem in the way i calculate the speed please can you give me a hint?
0 个评论
回答(1 个)
Torsten
2024-5-17
编辑:Torsten
2024-5-17
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!