咦？奇怪的知識又增加了！以數學建模的方式打開海嘯傳播模型（附部分matlab原始碼）

前言

之前關于森林火災蔓延模型小伙伴們反響都還不錯，今天我們對海嘯進行數學建模，

在文章初始，先喂自己袋鹽，

以下是博主精心整理的兩個matlab專欄，包含入門到精通及實戰內容，需要的小伙伴可根據自己需求自行訂閱，

MATLAB-30天帶你從入門到精通

https://blog.csdn.net/wenyusuran/category_10614422.html

MATLAB深入理解高級教程（附原始碼）

https://blog.csdn.net/wenyusuran/category_2239265.html

在博主的資源中也有各種演算法的應用實體源代碼，需要的小伙伴自取喲，

文中涉及到的代碼地址如下：

matlab原始碼集錦-海嘯傳播模型

https://download.csdn.net/download/wenyusuran/16237225

地震的形成是地殼在地球內動力地質作用下，地殼相鄰地塊發生緩慢的相對位移，地殼巖石發生變形，巖石應變數隨地塊逐步位移而增大，當巖石應變數大于巖石所能承受的強度時，巖石發生破裂，巖石彈性變形迅速恢復，這時儲存于巖石中的彈性變形能突然釋放，便形成地震，彈性變形地殼突然位移引起其上水體擾動，則形成涌浪、海嘯，所以，海嘯和地震在空間分布和發生時間上存在著密切的聯系，破壞型海嘯自然多分布在地震帶上，

我們現在都非常清楚海嘯是從深海開始的水波，通常是因為水下地震（雖然海嘯也可能是由水下滑坡或火山造成的），然后向海岸傳播，最初，海嘯的振幅相對較小（典型的是一米左右），這似乎使它們像風浪一樣無害，事實上，海嘯經常在深海中通過船只，甚至沒有人注意到，

海嘯傳播方式

由于海嘯波的傳播速度與海水深度的平方根有關，其傳播速度降低到每小時幾十公里，前進受到阻擋，就會形成十幾米和幾十米的浪高，沖向陸地.這種波長極長、速度極快的海嘯波，一旦從深海到達了岸邊，前進受到了阻擋，其全部的巨大能量，將變為巨大的破壞力量，摧毀一切可以摧毀的東西，造成巨大的災難，

海嘯是像地震一樣大的事件，海嘯的波長很大 - 典型的是200公里（與風浪相反，其波長通常接近100米），特別是，海嘯的波長遠遠大于海洋的深度（通常為2-3公里），因此，即使在深海中，海嘯的動態也基本上由淺水方程控制，這些方程的一個結果是海嘯的傳播速度，可以用公式近似：

是海洋的深度，是重力，因此，深水中的海嘯移動速度非常快 - 例如每小時500公里（每小時300英里）的速度是非常典型的; 例如，在不到一天的時間內就可以從日本旅行到美國，這是由于水的不可壓縮性（和質量守恒）; 非常寬和深的水波的巨大凈壓力（或者更準確地說，這種壓力的空間變化）迫使波浪的輪廓以很大的速度水平移動，（注意，這是海嘯波的相速度，而不是水分子本身的速度，它們要慢得多，）

當海嘯接近海岸時，深度當然會減少，導致海嘯減速，速度與深度的平方根成正比，如（1）所示，不幸的是，波浪淺灘然后迫使振幅以由格林定律控制的反向速率增加，

至少直到幅度變得與水深相當（此時，上述近似結果會被破壞;同樣，在兩個（水平）空間維度中，當海嘯向外擴散時，振幅會有一些衰減），如果假設海嘯開始，其初始振幅，處于深度，用比例關系計算振幅和深度，得海嘯的振幅（和水的深度），因此，例如，在2公里深度處初始幅度為1米的海嘯最終可以在岸邊附近達到約5米的最終幅度，同時仍然以每秒約10米的速度行進（每小時35公里，或22每小時英里數），我們現在就看到了當它到達岸邊時可能產生的影響，

雖然海嘯太大而無法控制（至少在深海），但我們至少可以在數學上對它們進行建模，從而可以高精度地預測它們在沿海各個地方的影響，用于執行此類模型的完整方程式和數值方法有些復雜，但是通過大量的簡化假設，就相對容易地得出一個已經預測了海嘯傳播基本特征的粗糙模型，如速度公式（1）和幅度比例法（2），我在下面給出了這個（標準）推導，這個論點在很大程度上是啟發式的; 實際上，下面的許多步驟都是非常有趣的分析問題，我們這里不做討論，

淺水波方程

當然，海洋是一種三維流體，但為了簡化分析，我們將考慮一個二維模型，其中唯一的空間變數是水平變數和垂直變數，其中是平衡海平面，我們用曲線模擬海底

因此海洋深度的測量位置，在任何時間和位置，水的高度（與海平面相比）將由未知的高度函式給出; 因此，在任何時候，海洋都占據了該地區

現在，我們通過在每個時間和海洋中的每個點分配速度矢量來模擬海洋內水的運動

我們做了不可壓縮性的基本假設，因此水的密度始終是恒定的，

根據牛頓第二定律，速度隨時間而變化，為了將此定律應用于流體，我們考慮沿著速度場流動的無窮小量的水，因此，在時間上，我們假設這個水量占據了我們所擁有的某些無限小的區域和某個位置

由于不可壓縮性，該區域保持不變，并且該無限小部分水的質量為，這個水體上會有兩股力量; 重力和壓力場的力，由壓力場給出，在海嘯的長度和時間尺度上，我們可以忽略其他力量的影響，如粘度或表面張力，然后牛頓定律給出了

這簡化了不可壓縮的歐拉方程

目前，沒有給出壓力，然而，我們可以通過假設（垂直）流體靜力平衡來簡化事物，即壓力的垂直效應抵消了重力的影響，我們還假設水面上的壓力為零，總之，這兩個假設迫使壓力成為靜水壓力

這反映了一個直觀合理的事實，即海洋下一點的壓力應該由該點以上的水的重量決定，

不可壓縮的歐拉方程現在簡化為

接下來我們將淺水近似表明水的波長遠大于水的深度，特別是，我們不期望變數中的速度場發生顯著變化，從而使ansatz成為可能

（這個ansatz應該用一粒鹽，特別是當應用于速度的分量時，實際上必須波動以適應海洋深度和高度函式的變化，但是，主要成分是速度是水平分量，在實際海嘯中，這確實表現為垂直不敏感，）

取的組分（4） ，和縮寫為，我們得到第一淺水波方程

這在泰勒擴展到第二個淺水波方程后簡化了

等式（6），（7）在未知數中是非線性的，然而，通過假設波的幅度與水的深度相比較小，可以近似線性化它們：

這個假設對于深海中的海嘯，甚至是中等深度的海嘯都是相當準確的，但是一旦海嘯到達海岸（動力學模型更難以模擬），當然是不合理的，

假設（11）已經簡化（7）到（近似）

至于（6），我們認為左側的第二個詞可以忽略不計，導致

為了解釋為什么是啟發式，我們希望這是情況，讓我們做的是擬設，并有幅度，分別是傳播相速度和波長; 讓我們也做出（合理的）假設，即空間變化比變化慢得多（即在波長范圍內大致恒定），因此我們可以（對于第一次近似）替換為，然后我們有

然后聯立等式（12）

從（11）我們期望，因此; 波的傳播速度比流體的速度快得多，特別是，我們預計會比小得多，這就解釋了為什么我們希望在（6）中第二個任期減少以獲得（13），

如果我們現在將上述ansatz插入（13），我們就得到了

將其與（14）相結合，我們已經得到了速度關系（1），

為了得到關系（2），我們必須更仔細地分析ansatz，首先，我們將（13）和（12）組合成高度函式的單個方程，實際上，在時間上區分（12）然后在（13）和（1）中代入

為了解決這個波動方程，我們使用標準的正弦曲線

是緩慢變化的功能，并且是一個小引數，插入此ansatz并提取項，我們得出了eikonal方程

和Hamilton-Jacobi方程

從eikonal方程我們看到它以速度傳播，假設向右傳播，我們就這樣做了

對于Hamilton-Jacobi方程，我們使用特征方法來解決它，乘以等式，我們得到

插入（15）和，獲得

這簡化為

因此，我們看到它在特征上是不變的，在另一方面，區分（15）中，我們看到

因此也是不變的特征，我們看到這是不變的特征，導致比例關系

給出（2），

原始碼

%Steven McHale
%Tsunami Model
%Shallow-Water Wave Equation
%Crank-Nicholson Discretization

clear;
clf;
clc;

% Constants
g  = 9.81;
u0 = 0;                               
v0 = 0;
b  = 0;                               
h0 = 5030;                             

% Define the x domain
ni = 151;                               
xmax = 100000;                      
dx = xmax/(ni-1);
x  = [0:dx:xmax];

% Define the y domain
nj = 151;                              
ymax = 100000;                      
dy = ymax/(nj-1);
y  = [0:dy:ymax];

% Define the wavespeed
wavespeed = u0 + sqrt(g*(h0 - b));

% Define time-domain
dt = 0.68*dx/wavespeed;
tmax = 1500;
%t = [0:dt:tdomain];
t=[1:dt:tmax];
courant = wavespeed*dt/dx;

% Build empty u, v, b matrices
u=zeros(length(x), length(y), length(t));
v=zeros(length(x), length(y), length(t));
b=zeros(length(x), length(y));

% Define h
h=zeros(length(x), length(y), length(t)); 
h(:,:,1) = 5000;                              
h((45000/100000*(length(x)-1)+1):floor(55000/100000*(length(x)-1)+1),(45000/100000*(length(y)-1)+1):floor(55000/100000*(length(y)-1)+1),1) = 5030;

%Define b
for i = 1:length(x)
    if x(i) > 20001
        b(:,i) = 0;
    elseif x(i) < 20000
        b(:,i) = 5000/20000*(20000-x(i));
    end
end

% Employ Lax
for n=1:(length(t)-1)
    for i=2:(ni-1)
        for j=2:(nj-1)    
            u(i,j,n+1) = ((u(i+1,j,n) + u(i-1,j,n) + u(i,j+1,n) + u(i,j-1,n))/4)...
                - 0.5*(dt/dx)*((u(i+1,j,n)^2)/2 - (u(i-1,j,n)^2)/2)...
                - 0.5*(dt/dy)*(v(i,j,n))*(u(i,j+1,n) - u(i,j-1,n)) - 0.5*g*(dt/dx)*(h(i+1,j,n)-h(i-1,j,n));
            
            v(i,j,n+1) = ((v(i+1,j,n) + v(i-1,j,n) + v(i,j+1,n) + v(i,j-1,n))/4)...
                - 0.5*(dt/dy)*((v(i,j+1,n)^2)/2 - (v(i,j+1,n)^2)/2)...
                - 0.5*(dt/dx)*(u(i,j,n))*(v(i+1,j,n) - v(i-1,j,n)) - 0.5*g*(dt/dy)*(h(i,j+1,n)-h(i,j-1,n));
            
            h(i,j,n+1) = ((h(i+1,j,n) + h(i-1,j,n) + h(i,j+1,n) + h(i,j-1,n))/4)...
                - 0.5*(dt/dx)*(u(i,j,n))*((h(i+1,j,n)-b(i+1,j)) - (h(i-1,j,n)-b(i-1,j)))...
                - 0.5*(dt/dy)*(v(i,j,n))*((h(i,j+1,n)-b(i,j+1)) - (h(i,j-1,n)-b(i,j-1)))...
                - 0.5*(dt/dx)*(h(i,j,n)-b(i,j))*(u(i+1,j,n)- u(i-1,j,n))...
                - 0.5*(dt/dy)*(h(i,j,n)-b(i,j))*(v(i,j+1,n) - v(i,j-1,n));

        end
    end

    % Define Boundary Conditions
    u(1,:,n+1) = 2.5*u(2,:,n+1) - 2*u(3,:,n+1) + 0.5*u(4,:,n+1);
    u(length(x),:,n+1) = 2.5*u(ni-1,:,n+1) - 2*u(ni-2,:,n+1) + 0.5*u(ni-3,:,n+1);
    u(:,1,n+1) = 2.5*u(:,2,n+1) - 2*u(:,3,n+1) + 0.5*u(:,4,n+1);
    u(:,length(y),n+1) = 2.5*u(:,nj-1,n+1) - 2*u(:,nj-2,n+1) + 0.5*u(:,nj-3,n+1);

    v(1,:,n+1) = 2.5*v(2,:,n+1) - 2*v(3,:,n+1) + 0.5*v(4,:,n+1);
    v(length(x),:,n+1) = 2.5*v(ni-1,:,n+1) - 2*v(ni-2,:,n+1) + 0.5*v(ni-3,:,n+1);
    v(:,1,n+1) = 2.5*v(:,2,n+1) - 2*v(:,3,n+1) + 0.5*v(:,4,n+1);
    v(:,length(y),n+1) = 2.5*v(:,nj-1,n+1) - 2*v(:,nj-2,n+1) + 0.5*v(:,nj-3,n+1);

    h(1,:,n+1) = 2.5*h(2,:,n+1) - 2*h(3,:,n+1) + 0.5*h(4,:,n+1);
    h(length(x),:,n+1) = 2.5*h(ni-1,:,n+1) - 2*h(ni-2,:,n+1) + 0.5*h(ni-3,:,n+1);
    h(:,1,n+1) = 2.5*h(:,2,n+1) - 2*h(:,3,n+1) + 0.5*h(:,4,n+1);
    h(:,length(y),n+1) = 2.5*h(:,nj-1,n+1) - 2*h(:,nj-2,n+1) + 0.5*h(:,nj-3,n+1);

end

figure(1)

%Animation of H wave propogating
for index=1:length(t)
    mesh(x,y,h(:,:,index))
    axis ([0 100000 0 100000 4990 5010])
    title ('AERSP 423 Computer Project Part II')
    xlabel('X Domain [m]')
    ylabel('Y Domain [m]')
    zlabel('Height [m]')
    pause(0.02)
end