필터 지우기
필터 지우기

Problem with my forward euler scheme

조회 수: 1 (최근 30일)
Vezzaz
Vezzaz 2022년 4월 28일
댓글: Vezzaz 2022년 4월 28일
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')

채택된 답변

Torsten
Torsten 2022년 4월 28일
편집: Torsten 2022년 4월 28일
%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')

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 MATLAB에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by