Effacer les filtres
Effacer les filtres

issue in calculation of wave speed

3 vues (au cours des 30 derniers jours)
am
am le 16 Mai 2024
Modifié(e) : am le 20 Mai 2024
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);
  1 commentaire
Walter Roberson
Walter Roberson le 17 Mai 2024
You have a problem with half_max_positions

Connectez-vous pour commenter.

Réponses (1)

Torsten
Torsten le 17 Mai 2024
Modifié(e) : Torsten le 17 Mai 2024
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);
  3 commentaires
Torsten
Torsten le 20 Mai 2024
Modifié(e) : Torsten le 20 Mai 2024
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
am le 20 Mai 2024
Modifié(e) : am le 20 Mai 2024
@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

Connectez-vous pour commenter.

Catégories

En savoir plus sur Mathematics dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by