Filter löschen
Filter löschen

issue in calculation of wave speed

26 Ansichten (letzte 30 Tage)
am
am am 16 Mai 2024 um 23:05
Bearbeitet: am am 20 Mai 2024 um 19:00
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 Kommentar
Walter Roberson
Walter Roberson am 17 Mai 2024 um 19:31
You have a problem with half_max_positions

Melden Sie sich an, um zu kommentieren.

Antworten (1)

Torsten
Torsten am 17 Mai 2024 um 19:36
Bearbeitet: Torsten am 17 Mai 2024 um 19:44
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 Kommentare
Torsten
Torsten am 20 Mai 2024 um 17:32
Bearbeitet: Torsten am 20 Mai 2024 um 18:07
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 am 20 Mai 2024 um 18:53
Bearbeitet: am am 20 Mai 2024 um 19:00
@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

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Mathematics finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by