HOW TO CREATE TSUNAMI MODEL
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
be discovered:
A = 5
X = 0 - 600
t = 0 - 60 minute
long = 0,5 - 30 kilometers
count the lamda, sigma and plot psi
4 commentaires
Sam Chak
le 17 Mai 2024
Déplacé(e) : Sam Chak
le 17 Mai 2024
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
Réponses (0)
Voir également
Catégories
En savoir plus sur Symbolic Computations in 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!