Azzera filtri
Azzera filtri

issue in calculation of wave speed

3 visualizzazioni (ultimi 30 giorni)
am
am il 16 Mag 2024
Modificato: am il 20 Mag 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);

Risposte (1)

Torsten
Torsten il 17 Mag 2024
Modificato: Torsten il 17 Mag 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 Commenti
Torsten
Torsten il 20 Mag 2024
Modificato: Torsten il 20 Mag 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 il 20 Mag 2024
Modificato: am il 20 Mag 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

Accedi per commentare.

Categorie

Scopri di più su Mathematics in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by