趣味で計算流砂水理

趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

数値計算とか河川工学とかプログラミングのことを書いています

MENU

粒度分布の関数形について

R

たまにはコアな土砂水理の話です。

本記事はGitHub、nbviewer、Colabでも公開しています。

※Colabの解説記事はこちら

a nbviewer Open In Colab

TL;DR(Bing Chat)
  • 粒度分布の関数形について：河床材料の粒度分布を対数正規分布型関数とTalbot型関数で近似し、フィッティングの結果を比較する。
  • 対象とする粒度分布：沖積河川上流域の礫床河川の交換層の粒度分布を用いる。
  • フィッティングの結果：対数正規分布型関数では粗めの成分と細かめの成分でフィッティングが異なる。Talbot型関数では指数を調整するだけで良くフィットする。

モチベーション

以下の2論文（ほぼ同じ内容）に記載される粒度分布の関数近似について、以前から気になっていたので少し考えてみました。

気になっていることは、粒度分布の関数形の一つとして論文に示されるTalbot型ですが、河床材料の場合、分級作用が働くため、平衡状態では必ずLog-normal型になるのではという点です。

少しだけ考えてみたのでまとめておきます。

本論文に関する予備知識

流砂水理の研究者なら知っていて当然ですが、一応整理しておきます。

本論文は、京都大学の藤田正治先生の研究チームがまとめた河床材料の粒度分布と空隙率の関係に関するシリーズの一部になります。 代表的な論文を整理しておきます。

この論文が最初の成果で粒径を球と仮定して数値実験により空隙率と粒度分布の関係をまとめています。

次が前述の論文で、実河川の粒度分布について近似関数を用いて分類し、関数形ごとの空隙率の特性をまとめています。

これまでに得られた粒度分布と空隙率の関係を用いて、一次元河床変動モデルを構築して数値実験を行っております。結論としてはなかなか難しいよねくらいの感じです。

一応、ここまでで研究は終わっているのですが、後に株式会社 建設技術研究所 岩見さんが発展させた研究を行っています。

全てはフォローできていませんが現時点では空隙率の変化を考慮した平面二次元河床変動計算を行っているようです。

以上が関連論文の整理になります。

著者の粒度分布の型の分類と疑問

論文には以下のように示されています。

Riverbed materials have a variety of different characteristic size of bed surface sediment (Bunte and Abt, 2001), but the grain size distribution could be classified into some types. One of the most typical density functions of grain size is a log-normal distribution. The grain size distribution of bed material in most sand-bed streams is unimodal and that in many gravel-bed rivers is bimodal (Parker, 2004). Also in mountainous rivers, the surface bed material has usually Talbot distribution of grain size (Tatsuzawa et al., 1998).

河床材料は、様々な特性を持つ粒径（Bunte and Abt, 2001）を持ちますが、粒径分布はいくつかのタイプに分類することができます。粒径の最も典型的な密度関数の一つは対数正規分布です。ほとんどの砂床河川の河床材料の粒径分布は単峰性で多くの礫床河川では双峰性(bimodal)です（Parker, 2004）。また、山地河川では、表面の河床材料は通常、粒径のタルボット分布を持っています（Tatsuzawa et al., 1998）。

のなかで「山地河川では、表面の河床材料は通常、粒径のタルボット分布を持っています。」部分が引っかかっています。 山地河川でも分級作用があるので、対数正規分布になるのでは？と思っています。

参照文献の

竜澤 宏昌, 林 日出喜, 長谷川 和義 : 渓流河川における河床砂礫の混合特性と階段状河床形の形状特性, 水工学論文集42巻, 1998

をみると、渓流河川の粒度分布はTalbot型関数で表現できるとありますが、分級機構については十分に検討されていません。

実河川の粒度分布での検証

対象とする粒度分布

ある沖積河川上流域の比較的急流な礫床河川の粒度分布を対象とします。 表層を剥いだ下層、いわゆる交換層の粒度分布になります。

粒度分布をみると最大粒径が100mmとなっていますが、これは材料採取の際に100mm以上を採取していないため?でしょうか。

この粒度分布を対象に各種関数でフィッティングしていきます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, lognorm

d = np.array([100,75,53,37.5,26.5,19,9.5,4.75,2,0.85,0.425,0.25,0.106,0.075,0.053])
p = np.array([100,94.6,77.9,68,55.9,49.2,37.2,30.1,23.1,14.3,6.4,2.8,0.7,0.5,0])
d = d[::-1]
p = p[::-1]

%matplotlib inline
fig, ax = plt.subplots()

ax.plot(d, p, marker='o')
ax.set_xscale('log')
ax.grid(which='major',color='black',linestyle='-')
ax.grid(which='minor',color='gray',linestyle='--')
ax.set_xlabel('diameter[mm]')
ax.set_ylabel('percentage of grain size[%]')

plt.show()

対数正規分布型関数によるフィッティング

対数正規分布の確率密度関数は以下に定義される。

\begin{align} f(x) = \frac{1}{\sqrt{2\pi}\sigma x} \exp\left(-\frac{(\log x - \mu)^2}{2\sigma^2}\right) \end{align}

ここに、\mu:期待値、\sigma：標準偏差とする。

また、対数正規分布の累積分布関数は以下に定義される。

\begin{align} F(x) = \frac{1}{2} \operatorname{erfc} \left(-\frac{\ln x - \mu}{\sqrt{2}\sigma}\right) \end{align}

ここに、\operatorname{erfc}(x)は相補誤差関数を表す。

これらより、\mu、\sigmaの2パラメータでフィッティングする。

面倒なので、scipy.stats.lognormを使います。

引数の定義が正規分布の場合と異なるのでわからない方は以下を参照下さい。

フィッティング結果は次図のとおりです。

粒度分布が複雑な形状なので、粗めの成分でフィッティングする場合と細かめの成分でフィッティングする場合の2ケースが考えられます。 言い換えると対数正規分布でのフィッティングが難しいと考えられます。 bimodal型に近いかもしれませんね。フィッティングする場合は、対数正規分布の重ね合わせでできそうです。後日やってみます。

%matplotlib inline

fig, ax = plt.subplots()

ax.plot(d, p, marker='o', label='observed')

xd = np.arange(-2, 3.5, 0.1)
y = lognorm.cdf(10**xd, s=1, scale=np.exp(3.2), loc=0)
ax.plot(10**xd, y*100, label="$\sigma=1, \mu=24.53$")

y = lognorm.cdf(10**xd, s=2.5, scale=np.exp(3), loc=0)
ax.plot(10**xd, y*100, label="$\sigma=2.5, \mu=20.09$")

ax.set_xscale('log')
ax.grid(which='major',color='black',linestyle='-')
ax.grid(which='minor',color='gray',linestyle='--')
ax.set_xlabel('diameter[mm]')
ax.set_ylabel('percentage of grain size[%]')
ax.legend()

plt.show()

Talbot型関数によるフィッティング

そもそものTalbot型の粒度分布について簡単にまとめておきます。

フラー(W.B.Fuller)とトンプソン(S.E.Thompson)によって確立された(1907)、粒子を最密充填にするための粒度分布曲線(Fuller曲線)が次式で示されます。 最大密度を与える理想粒度としてコンクリート骨材分野で提案された式らしいです。

\begin{align} f(x) = \sqrt{\frac{x}{d_{max}}} \end{align}

上式の一般形が下式のTalbot式になります。(1923)

\begin{align} f(x) = \left(\frac{x}{d_{max}}\right)^n \end{align}

著者は河床材料にフィッティングしやすいようにTalbot式を以下のように修正しています。

\begin{align} f(x) = \left(\frac{\log{x}-\log{d_{min}}}{\log{d_{max}}-\log{d_{min}}}\right)^{n_T} : n_T>1 \end{align}

この式を使ってフィッティングしていきます。パラメータはn_Tのみです。

参考文献

福本武明, 土の連続型粒度式, 土木学会論文集1994

フィッティング結果は次図のとおりです。

指数を調整するだけでは、良くフィッティングできています。 これだけフィットすること考えると、Talbot型の粒度分布が形成される物理的要因がありそうです。

%matplotlib inline

n = 0.45
f = (d/d[-1])**n

nt = 2.7
ft = ((np.log10(d)-np.log10(d[0]))/(np.log10(d[-1])-np.log10(d[0])))**nt

fig, ax = plt.subplots()
ax.plot(d, p, marker='o', label='observed')
ax.plot(d, f*100, label='Talbot $n=0.45$')
ax.plot(d, ft*100, label='Talbot_modified $n\_T=2.7$')
ax.set_xscale('log')
ax.grid(which='major',color='black',linestyle='-')
ax.grid(which='minor',color='gray',linestyle='--')
ax.set_xlabel('diameter[mm]')
ax.set_ylabel('percentage of grain size[%]')
ax.legend()
plt.show()

まとめ

  • 実河川の河床材料をみても、Talbot型の粒度分布は存在しそう。
  • Talbot型の粒度分布が形成される要因については不明である。出水時には常に全粒径が動くとか？
  • 数値実験をやって検証しようと思います。

参考図書

勉強会の議事メモ

R

勉強会お疲れ様でした。

粒度分布の関数形について

Talbot型の物理的意味について整理

computational-sediment-hyd.hatenablog.jp

単位流域から流出する流量波形のガンマ関数近似について

理論の整理

Mann-Kendall検定の問題点

理論の整理、問題点の整理

一般座標系平面二次元河床変動計算を書く

久しぶりに書きましょう。我々のモデルが無いのは不便

久しぶりに本質的な議論ができて楽しかったです。

実務では十分にできない物事を細部まで論理的に分析することは重要ですね。

conda install時にOpenSSLのエラーが出たときの解決方法(windows10環境)

R

condaの環境を整理した後にconda installを実行した時に以下のOpenSSLに関するエラーが発生したため、その解決方法です。 情報があまりなかったので簡単にまとめておきます。

Collecting package metadata (current_repodata.json): failed

CondaSSLError: OpenSSL appears to be unavailable on this machine. OpenSSL is required to
download and install packages.

Exception: HTTPSConnectionPool(host='conda.anaconda.org', port=443): Max retries exceeded with url: /conda-forge/win-64/current_repodata.json (Caused by SSLError("Can't connect to HTTPS URL because the SSL module is not available."))

参照サイト

community.anaconda.cloud

github.com

参照サイトによると、OpenSSL関連のライブラリが上手く参照できていないことが原因のようです。 そのため、シンボリックリンクを使って以下のコマンドを実行することで解決できます。もちろんコピーしても大丈夫です。

mklink "C:\ProgramData\Anaconda3\DLLs\libcrypto-1_1-x64.dll" "C:\ProgramData\Anaconda3\Library\bin\libcrypto-1_1-x64.dll"
mklink "C:\ProgramData\Anaconda3\DLLs\libssl-1_1-x64.dll" "C:\ProgramData\Anaconda3\Library\bin\libssl-1_1-x64.dll"

pythonモジュール：rioxarrayをアップデートしたら結構変わっていたのでよく使うところをまとめました

R

このブログではおなじみのpythonモジュール：rioxarrayを最新版にアップデートすると、主要なメソッドが変わっていたのでまとめておきます。

rioxarrayを使っている主な記事は以下です。

GeoTIFFファイルのreadが廃止：xarray.open_rasterio ⇒ rioxarray.open_rasterio

stackoverflow.com

xarray.open_rasterioは、xarrayの2023.04の更新で完全に廃止されました。

open_rasterioはrioxarrayモジュールが無いと動かないので、わかりやすくなったかと思います。

以前は、

import xarray as xr
ds = xr.open_rasterio('A.tif')

でも、readできましたが、今後は、

import rioxarray
ds = rioxarray.open_rasterio('A.tif')

のみになります。

crsの確認方法

github.com

上のページのとおりですが、以前は、

cog_url = (
    "https://oin-hotosm.s3.amazonaws.com/"
    "5d7dad0becaf880008a9bc88/0/5d7dad0becaf880008a9bc89.tif"
)
import xarray as xr
xr.open_rasterio(cog_url).crs

# '+init=epsg:3857'

で座標情報を確認できましたが、rioxarrayを用いた場合は、

import rioxarray
rioxarray.open_rasterio(cog_url).spatial_ref

# <xarray.DataArray 'spatial_ref' ()>
# array(0)
# Coordinates:
#     spatial_ref  int64 0
# Attributes:
#     spatial_ref:  PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WG...
#     crs_wkt:      PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WG...

のように、spatial_refというCoordinatesを作成されて、そのAttributesに座標情報が格納されています。

epsgを確認する場合は、

ds = rioxarray.open_rasterio(cog_url)

ds.rio.crs.to_proj4()
# '+init=epsg:3857'

ds.rio.crs.to_epsg()
# 3857

などになります。

数値解析の基礎：10分解説

R

個人的な勉強会で作った資料です。 詳しい方から見たらふざけるな的な内容ですが10分解説なので許して下さい。

本記事はGitHub、nbviewer、Colab、Marp（スライド形式）でも公開しています。

※Colabの解説記事はこちら

※Marpの解説記事はこちら

a nbviewer Open In Colab a

数値解析のアプローチ

  • 自然現象や社会現象など世の中の多くの現象は微分方程式によって表すことができる（当然その限りではない）。
  • 数値解析とは微分方程式によってモデル化された現象をコンピュータで取り扱える形にして計算を行なうことを示す。
  • 微分方程式を離散化することにより代数方程式を導出する。離散化手法のうち 、流体解析でよく使われる差分法について説明する

※ 茨城大学 田中伸厚先生の講義資料を引用 : http://www.mech.ibaraki.ac.jp/~tanaka/Lecture/sim1/part1.pdf

離散化について -差分法を例に-

ある現象が次の常微分方程式で表すことができるとする。

\dfrac{dN(t)}{dt} = -k N(t)

ここに、k:崩壊定数とする。

上式は連続系でありそのままコンピュータで取り扱うことができない。

そこで上式を

\dfrac{N(t+\Delta t)-N(t)}{\Delta t} = -k N(t)

のようにコンピュータ上で取り扱える式（代数方程式＝四則演算のみの式）に変換することを離散化（差分法）という。

また、代数方程式は添字を用いて以下のように記述することが多い。以降は以下の表記を使用する。

\dfrac{N^{n+1}-N^n}{\Delta t} = -k N^n

計算例1：放射性物質の崩壊

現象

放射性物質の原子核数の時間変化

出典：https://studyphys.com/half-life/

数学モデル（微分方程式）

以下の常微分方程式で表すことができる。

\dfrac{dN(t)}{dt} = -k N(t)

ここに、k:崩壊定数とする。

数値モデル（代数方程式）

差分法を用いると以下のとおりとなる。

\begin{align} \dfrac{N^{n+1}-N^n}{\Delta t} &= -k N^n \\ N^{n+1} &= N^n - \Delta t k N^n \end{align}

この式よりN^0の値（初期条件）を与えて逐次計算を行うことにより数値解を得ることができる。

厳密解

前出の常微分方程式は厳密解（数値計算を用いずに解析的に得られる解。解析解とも呼ばれる）を求めることができる。

導入過程は省略するとが、変数分離を用いて式変形することにより下式が得られる。

\begin{align} N(t) = N^0 e^{-kt} \end{align}

本式による計算結果と数値計算の結果を比較することにより計算の妥当性を評価することができる。

計算結果

導出した代数方程式を用いて数値解析を実施する。 \Delta tを0.1,0.5,1.0,1.5と変化させた場合の計算結果と厳密解との比較を示す。 計算条件は、N^0は100、kは1とする。

import numpy as np
import matplotlib.pyplot as plt

N0 = 100
k = 1.0
exact = lambda T: N0 * np.exp(-k*T)

def simulation(T):
    N = np.zeros_like(T)
    N[0] = N0
    for i in range(1, len(N)):
        N[i] = N[i-1] - dt*k*N[i-1] 
    return N

def mkfig(ax, T):
    ax.plot(T, exact(T),label='exact')
    ax.plot(T, simulation(T),label='simulation')
    ax.set_xlabel('T')
    ax.set_ylabel('N')
    ax.legend(loc = 'upper right')

\Delta t = 0.1の計算結果

dt = 0.1
T = np.arange(0, 20, dt)
fig, ax = plt.subplots()
mkfig(ax, T)

png

\Delta t = 0.5の計算結果

dt = 0.5
T = np.arange(0, 20, dt)
fig, ax = plt.subplots()
mkfig(ax, T)

png

\Delta t = 1.0の計算結果

dt = 1.0
T = np.arange(0, 20, dt)
fig, ax = plt.subplots()
mkfig(ax, T)

png

\Delta t = 1.5の計算結果

dt = 1.5
T = np.arange(0, 20, dt)
fig, ax = plt.subplots()
mkfig(ax, T)

png

  • \Delta tによって解析結果が異なる。
  • \Delta t=0.1は、厳密解と数値解析による差はほとんどないが、\Delta tが大きくなるにつれて差が大きくなる。\Delta t=1.5は異常な結果を示している。
  • 差分法では\Delta tを小さく設定することが望ましいが、計算資源の制約で限界がある。

他の差分法

  • 今回はもっと基本的な前進差分スキーム（nの値からからn+1の値を求める）を使用したが、後退差分スキーム（n+1の値からからn+1の値を求める）も存在する。
  • 時間発展の前進差分スキームは陽解法、後退差分スキームは陰解法と呼ばれる。
\begin{align} \dfrac{N^{n+1}-N^n}{\Delta t} &= -k N^{n+1} \\ N^{n+1} &= \dfrac{N^n}{1+ \Delta t k} \end{align}
  • 以下に示す計算結果のとおり、\Delta t=1.5でもそれなりの数値計算結果が得られることがわかる。
  • 陰解法はメリットしか無いように思われるが、陰解法では上記のような単純な式形になることは稀であり、ほとんどのケースで複雑な式形になり、数値解を得るために収束計算が必要となる。
import numpy as np
import matplotlib.pyplot as plt

N0 = 100
k = 1.0
exact = lambda T: N0 * np.exp(-k*T)

def simulation2(T):
    N = np.zeros_like(T)
    N[0] = N0
    for i in range(1, len(N)):
        N[i] = N[i-1]/(1+dt*k)
    return N

def mkfig2(ax, T):
    ax.plot(T, exact(T),label='exact')
    ax.plot(T, simulation2(T),label='simulation')
    ax.set_xlabel('T')
    ax.set_ylabel('N')
    ax.legend(loc = 'upper right')

\Delta t = 1.5の計算結果

dt = 1.5
T = np.arange(0, 20, dt)
fig, ax = plt.subplots()
mkfig2(ax, T)

png

計算例2：移流方程式

現象

移流とは下図のようにある物理量が流れによって運ばれる現象を示す。この現象を計算対象とする。

import numpy as np
import matplotlib.pyplot as plt

c = 1.0
X = np.linspace(-5,5,101)
exact = lambda t: np.exp(-(X-c*t)**2)
fig, ax = plt.subplots()

ax.plot(X, exact(0),label='0 sec')
ax.plot(X, exact(2),label='2 sec')
ax.set_xlabel('x')
ax.set_ylabel('u')
ax.legend(loc = 'upper right')

png

数学モデル（微分方程式）

以下の偏微分方程式で表すことができる。

\begin{align} \dfrac{\partial u(t,x)}{\partial t} + c\dfrac{\partial u(t,x)}{\partial x} =0 \end{align}

数値モデル（代数方程式）

今回は3つの差分スキームを使用する。

前進差分

\begin{align} \dfrac{u^{n+1}_i-u^n_i}{\Delta t} &= -c\dfrac{u^n_{i+1}-u^n_{i}}{\Delta x} \\ u^{n+1}_i &= u^n_i - \dfrac{c\Delta t}{\Delta x}\ (u^n_{i+1}-u^n_{i}) \end{align}

後退差分

\begin{align} \dfrac{u^{n+1}_i-u^n_i}{\Delta t} &= -c\dfrac{u^n_{i}-u^n_{i-1}}{\Delta x} \\ u^{n+1}_i &= u^n_i - \dfrac{c\Delta t}{\Delta x}\ (u^n_{i}-u^n_{i-1}) \end{align}

中心差分

\begin{align} \dfrac{u^{n+1}_i-u^n_i}{\Delta t} &= -c\dfrac{u^n_{i+1}-u^n_{i-1}}{2\Delta x} \\ u^{n+1}_i &= u^n_i - \dfrac{c\Delta t}{\Delta x}\dfrac{u^n_{i+1}-u^n_{i-1}}{2} \end{align}

この式よりu^0の値（初期条件,空間分布を持つ）を与えて逐次計算を行うことにより数値解を得ることができる。 なお、境界条件も必要であるが今回は境界付近を計算しないため0とする。

厳密解

一次元の移流方程式は条件によって厳密解を得ることができる。 u(0,x)=e^{-x^2}を与えると厳密解として次式が得られる。

\begin{align} u(t,x) = e^{-(x-ct)^2} \end{align}

計算結果

導出した代数方程式を用いて数値解析を実施する。 計算条件は、c=1、\Delta t=0.1、\Delta t=0.01として、厳密解との比較を行なう。

import numpy as np
import matplotlib.pyplot as plt

c = 1.0
X = np.linspace(-5,5,101)
dx = 0.1
exact = lambda t: np.exp(-(X-c*t)**2)

def s01(dt): #前進差分
    unew = exact(0)
    for _ in range(int(2/dt)):
        u = unew.copy()
        for i in range(1,len(u)-1):
            unew[i] = u[i] - c * dt/dx * (u[i+1]-u[i])
    return unew

def s02(dt): #後退差分
    unew = exact(0)
    for _ in range(int(2/dt)):
        u = unew.copy()
        for i in range(1,len(u)-1):
            unew[i] = u[i] - c * dt/dx * (u[i]-u[i-1])
    return unew

def s03(dt): #中心差分
    unew = exact(0)
    for _ in range(int(2/dt)):
        u = unew.copy()
        for i in range(1,len(u)-1):
            unew[i] = u[i] - c * dt/dx * (u[i+1]-u[i-1])/2
    return unew

前進差分

適切な解を得ることができない。

fig, ax = plt.subplots()
ax.plot(X, exact(0),label='0 sec')
ax.plot(X, exact(2),label='exact:2sec')
ax.plot(X, s01(0.01),label='simulation:2sec dt=0.01')
ax.set_xlabel('x')
ax.set_ylabel('u')
ax.legend(loc = 'best')

png

後退差分

波形が拡散してピークが減衰する。

fig, ax = plt.subplots()
ax.plot(X, exact(0),label='0 sec')
ax.plot(X, exact(2),label='exact:2sec')
ax.plot(X, s02(0.01),label='simulation:2sec dt=0.01')
ax.set_xlabel('x')
ax.set_ylabel('u')
ax.legend(loc = 'best')

png

中心差分

若干の差異はあるが他の差分スキームと比較すると厳密解に近い結果になっている。

fig, ax = plt.subplots()
ax.plot(X, exact(0),label='0 sec')
ax.plot(X, exact(2),label='exact:2sec')
ax.plot(X, s03(0.01),label='simulation:2sec dt=0.01')
ax.set_xlabel('x')
ax.set_ylabel('u')
ax.legend(loc = 'best')

png

  • 差分スキームによって解析結果が異なる。
  • 前進差分は\Delta tをいくら小さくしても適切な解を得ることができない。
  • 後退差分は拡散するため、\Delta t, \Delta xを十分に小さくとる必要がある。
  • 今回のケースでは中心差分が良さそうであるが、適切な解が得られない場合が多い。（発散し易い。）
  • その場合、より高精度なスキームが必要となり代数方程式の導出が複雑となる。⇒ 不定流計算などでは煩雑な代数方程式が必要

陰解法

一次元移流方程式の陰解法は前出の常微分方程式と比べてかなり複雑な式形となる。

例として陰解法で中心差分の代数方程式を導出すると次式となる。

\begin{align} \dfrac{u^{n+1}_i-u^n_i}{\Delta t} &= -c\dfrac{u^{n+1}_{i+1}-u^{n+1}_{i-1}}{2\Delta x} \\ u^{n+1}_i + \dfrac{c\Delta t}{\Delta x}\dfrac{u^{n+1}_{i+1}-u^{n+1}_{i-1}}{2} &= u^n_i \end{align}

上式のみでは、n+1の項の値を求めることができない。そのため、他の空間離散点（i,i+1,i+2など）の式を全て連立させた次式より解を計算する。

\begin{pmatrix} \ddots & \ddots & & & \\ \dfrac{c\Delta t}{2\Delta x} & 1 & -\dfrac{c\Delta t}{2\Delta x} & & \huge{0} \\ & \dfrac{c\Delta t}{2\Delta x} & 1 & -\dfrac{c\Delta t}{2\Delta x} & \\ \huge{0} & & \dfrac{c\Delta t}{2\Delta x} & 1 & -\dfrac{c\Delta t}{2\Delta x} \\ & & & \ddots & \ddots \end{pmatrix} \begin{pmatrix} \vdots \\ u^{n+1}_{i+1} \\ u^{n+1}_i \\ u^{n+1}_{i-1} \\ \vdots \end{pmatrix} = \begin{pmatrix} \vdots \\ u^{n}_{i+1} \\ u^{n}_i \\ u^{n}_{i-1} \\ \vdots \end{pmatrix}

左辺の対角行列の逆行列を計算することにより、n+1の項を求めることができる。

ここでは計算は省力するが、陰解法は陽解法と比べて煩雑な計算が必要となることが理解できる。

補足として、上式の左辺の1つ目の行列は三重対角行列となっており、TDMA法を用いて少ない計算ステップ数で逆行列を求めることができる。この方法は不等流計算の平均流速公式レベル3の解法でも使用する。

離散化方法のポイント

一般的に精度と安定性は両立しない。目的に応じて差分スキームを選択する必要がある。

精度

  • 離散化による打ち切り誤差の大きさ(○次精度のスキームなどと呼ぶ)

安定性

  • 計算誤差が時間の経過とともに発展しないこと。
    • 中心差分は計算誤差が発展しやすく不安定なスキームである。
    • 陰解法は安定性が高い。

まとめ

  • 数値解析は微分方程式を差分法などにより離散化した代数方程式から数値解を得ることである。
  • 差分法では、\Delta x, \Delta tの大きさによって解析結果が大きく異なるため、適切な条件設定が必要である。
  • 差分法には様々な方法（スキーム）があり、目的に応じた選択が必要である。

⇒ 今回の勉強で実施する不等流計算、貯留関数法では、差分スキームはそれほど重要ではないが\Deltaの大きさについてはある程度考慮する必要がある。

⇒ 不定流計算などを実施する場合は、差分スキームが重要

参考資料

備忘録：windowsでHEICをJPGに変換する。

R

iphoneでおなじみのHEICファイルですが、他ソフトでの使用のためJPGファイルに変化する機会は多いです。

いつもはGIMPなどを使っていましたが、 Windowsでは、Microsoft Storeの「Image Converter」で一括変換できることを知りました。

参照：【2023年最新】HEICファイルをJPGに変換する方法