issue in calculation of wave speed
3 次查看(过去 30 天)
显示 更早的评论
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
% 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 个)
Torsten
2024-5-17
编辑:Torsten
2024-5-17
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 个评论
Torsten
2024-5-20
编辑:Torsten
2024-5-20
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
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!