Effacer les filtres
Effacer les filtres

Problem with my forward euler scheme

1 vue (au cours des 30 derniers jours)
Vezzaz
Vezzaz le 28 Avr 2022
Commenté : Vezzaz le 28 Avr 2022
I am trying to make a code to run the explicit forward euler scheme to solve the 1d wave equation(Ut=c*Ux) with no boundary conditions and I cant get the indexing right and it is leading to an error. The last term where it is j-1 is the problem. It wont run unless I switch it to j= 2:Nt but that is leaving me with an incorrect answer. If i left it as j the answer still doesnt seem right but it at least looks like a plot that would be close. It might not be the indexing and my code could just be completely wrong, but I have no idea. Any help would be greatly appreciated.
clear all
close all
%parameters
Nx = 50; % Number of spatial steps
xl = 0; % Left boundary
xr = 4; % Right boundary
dx = 0.04; % Spatial step
x = xl:dx:xr; % x
dt=0.02; % time step
tf=2; % final time
Nt = round(tf/dt); % Number of time steps
t = 0:dt:2; % t
c=0.5; %CFL condition
u=zeros(length(x)); %initialize u
%initial condition
for i = 1:Nx+1
if x<0
u(i,1) = 0; % u(x,0)
else
u(i,1) = 1; % u(x,0)
end
end
%explicit forward euler
for j = 1:Nt % Time Loop
for i = 2:Nx
u(i+1,j)=u(i,j)+c*(dt/(2*dx))*(u(i,j+1)-u(i,j-1));
end
end
plot(u(:,1))
xlabel('Time')
ylabel('Solution')
title('Forward Euler Scheme, part a, dx=0.04 dt=0.02')

Réponse acceptée

Torsten
Torsten le 28 Avr 2022
Modifié(e) : Torsten le 28 Avr 2022
%parameters
c = 0.05; %Velocity
Nx = 50; % Number of spatial steps
xl = 0; % Left boundary
xr = 4; % Right boundary
x = linspace(xl,xr,Nx) ; % x
dx = (xr-xl)/(Nx-1); % Spatial Step
dt = 0.9*dx/c; % CFL condition
tf=20; % final time
t = 0:dt:tf; % t
Nt = numel(t); % Number of time steps
u=zeros(numel(x),numel(t)); %initialize u
%initial condition
u(2:nx,1) = 1.0;
%explicit forward euler
for j = 1:Nt-1 % Time Loop
u(2:Nx,j+1) = u(2:Nx,j)-c*dt/dx*(u(2:Nx,j)-u(1:Nx-1,j));
end
plot(x,[u(:,10),u(:,end)])
xlabel('Time')
ylabel('Solution')
title('Forward Euler Scheme, part a, dx=0.04 dt=0.02')

Plus de réponses (0)

Catégories

En savoir plus sur MATLAB 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