フィルターのクリア

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')
  1 件のコメント
Vezzaz
Vezzaz 2022 年 4 月 28 日
Thank you!

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeMATLAB についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by