Clear Filters
Clear Filters

Problem with my forward euler scheme

1 view (last 30 days)
Vezzaz
Vezzaz on 28 Apr 2022
Commented: Vezzaz on 28 Apr 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')

Accepted Answer

Torsten
Torsten on 28 Apr 2022
Edited: Torsten on 28 Apr 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')

More Answers (0)

Categories

Find more on MATLAB in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!