HOW TO CREATE TSUNAMI MODEL

5 次查看（过去 30 天）
meoui
meoui 2024-5-17
移动：Sam Chak 2024-5-17
be discovered:
A = 5
X = 0 - 600
t = 0 - 60 minute
long = 0,5 - 30 kilometers
count the lamda, sigma and plot psi
  4 个评论
显示 2更早的评论
Sam Chak
Sam Chak 2024-5-17
移动：Sam Chak 2024-5-17
Ran in:
Hi @Fadilla Neisa
Looking for this?
%% Tsunami model
syms L H depthratio g positive
syms x t w T R U(x)
L1 = depthratio*L;
L2 = L;
h1 = depthratio*H;
h2 = H;
h(x) = x*H/L;
c1 = sqrt(g*h1);
c2 = sqrt(g*h2);
u(x,t) = U(x)*exp(1i*w*t);
u1(x,t) = T*exp(1i*w*(t + x/c1));
u2(x,t) = exp(1i*w*(t + x/c2)) + R*exp(1i*w*(t - x/c2));
wavePDE(x,t) = diff(u, t, t) - g*diff(h(x)*diff(u, x), x);
slopeODE(x) = wavePDE(x, 0);
U(x) = dsolve(slopeODE);
Const = setdiff(symvar(U), sym([depthratio, g, H, L, x, w]));
du1(x) = diff(u1(x, 0), x);
du2(x) = diff(u2(x, 0), x);
dU(x) = diff(U(x), x);
eqs = [ U(L1) == u1(L1, 0), U(L2) == u2(L2, 0),...
dU(L1) == du1(L1), dU(L2) == du2(L2)];
unknowns = [Const(1), Const(2), R, T];
[Cvalue1, Cvalue2, R, T] = solve(eqs, unknowns);
U(x) = subs(U(x), {Const(1), Const(2)}, {Cvalue1, Cvalue2});
%% Parameters
g = 9.81;
L = 2;
H = 1;
depthratio = 0.04;
h1 = depthratio*H;
h2 = H;
L1 = depthratio*L;
L2 = L;
c1 = sqrt(g*h1);
c2 = sqrt(g*h2);
A = 0.3;
%% incoming soliton wave
soliton = @(x,t) A.*sech(sqrt(3/4*g*A/H)*(x/c2+t)).^2;
%% creating time scale and sample points
Nt = 64;
TimeScale = 40*sqrt(4/3*H/A/g);
W = [0, 1:Nt/2 - 1, -(Nt/2 - 1):-1]'*2*pi/TimeScale;
Nx = 100;
x_min = L1 - 4;
x_max = L2 + 12;
X12 = linspace(L1, L2, Nx); % slope region
X1 = linspace(x_min, L1, Nx); % shallow water region
X2 = linspace(L2, x_max, Nx); % deep water region
%% Fourier transform of the incoming soliton
S = fft(soliton(- 0.8*TimeScale*c2, linspace(0, TimeScale, 2*(Nt/2) - 1)))';
S = repmat(S, 1, Nx);
%% Construct a traveling wave solution based on S
soliton = real(ifft(S.*exp(1i*W*X2/c2)));
%% Construct a reflected wave
R = double(subs(subs(vpa(subs(R)), w, W), x ,X2));
R(1,:) = double((1 - sqrt(depthratio)) / (1 + sqrt(depthratio)));
reflectedWave = real(ifft(S.*R.*exp(-1i*W*X2/c2)));
%% Construct a transmitted wave
T = double(subs(subs(vpa(subs(T)), w, W), x, X1));
T(1,:) = double(2/(1+sqrt(depthratio)));
transmittedWave = real(ifft(S.*T.*exp(1i*W*X1/c1)));
%% Construct the wave at the slope region
U12 = double(subs(subs(vpa(subs(U(x))), w, W), x, X12));
U12(1,:) = double(2/(1 + sqrt(depthratio)));
U12 = real(ifft(S.*U12));
%% Plot the Solution
soliton = interpft(soliton, 10*Nt);
reflectedWave = interpft(reflectedWave, 10*Nt);
U12 = interpft(U12, 10*Nt);
transmittedWave = interpft(transmittedWave, 10*Nt);
figure('Visible', 'on');
plot([x_min, L1, L2, x_max], [-h1, -h1, -h2, -h2], 'linewidth', 2, 'Color', '#BEA256')
axis([x_min, x_max, -H-0.1, 0.6]), grid on
hold on
line1 = plot(X1, transmittedWave(1,:), 'linewidth', 4, 'Color', '#8CD0E1');
line12 = plot(X12, U12(1,:), 'linewidth', 3, 'Color', '#2AA7C0');
line2 = plot(X2, soliton(1,:) + reflectedWave(1,:), 'linewidth', 2, 'Color', '#0A7B88');
text(6, -0.6, 'Deep Water region')
for t = 2:size(soliton, 1)*0.35
line1.YData = transmittedWave(t,:);
line12.YData = U12(t,:);
line2.YData = soliton(t,:) + reflectedWave(t,:);
pause(0.1)
end

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Applications 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by