​
Slide
1
of
18

微分方程式を解こう​〜 Mathematica 12での

Wolfram Research 丸山耕司
kojim@wolfram.com
​​​​
ウェビナーではこのノートブックに式を追加したり変更したりしながら説明いたします．そのため，記述が足りなかったり，Shift + Enterしただけでは正しく動作しない部分も一部あることをご承知おきください．

​

​
Slide
2
of
18

今日のスピーカー : 丸山耕司

◼
  • Sales Engineer at Wolfram Research
    客員教員 at 大阪市立大学，大阪大学
    非常勤講師 at 早稲田大学，芝浦工業大学
    • ◼
  • Articles in Wolfram Blog : https://blog.wolfram.com/author/koji-maruyama/
    ・機械学習機能を利用した心不全による死亡のリスク予測
    ・Distinguishing Risks of Modes of Cardiac Death in Heart Failure with Machine Learning
    ・非線形偏微分方程式への有限要素法の適用
    • ◼
  • 専門: 理論物理(量子情報，量子制御理論，etc.)
    • RadFan : 放射線医学に関する業界向け情報誌．2021年6月号に「Wolfram言語/
    (R)
    Mathematica
    による画像処理・解析 」を寄稿．
    ​
    Slide
    3
    of
    18

    DSolve / DSolveValueで微分方程式を解く

    DSolveの使い方 (DSolveValueも同様) :DSolve[ 式(等号は==のようにふたつ並べる！) , 従属変数 , パラメータ ]
    DSolve[{eqn,cond},y[x],x]​​
    式，初期条件を変数にまるごと代入
    グラフ等を描くときのために, DSolveの結果をまるごとなにかの変数に代入してしまうと便利
    solution=DSolve[{y''[t]+1/10y'[t]+y[t]0,y[0]1,y'[0]0},y[t],t]
    解がy→...の形になっている. これは, /.solutionとすると，直前のyを矢印以下の表現で置き換えて評価せよ, という意味になる.
    Plot[y[t],{t,0,20}]
    DSolveValueを使うと結果の表式だけを出力として得られるので，この方が単純でわかりやすいことも多い．
    以下，DSolveとDSolveValueは言葉としては区別しない．
    solution2=DSolveValue[{y''[t]+1/10y'[t]+y[t]0,y[0]1,y'[0]0},y[t],t]
    Plot[solution2,{t,0,20}]
    例題1a:微分方程式y'=
    2
    x
    yを解いてみよう
    例題1b:同じy'=
    2
    x
    yを初期条件y(0)=2の下で解いてみよう
    例題1c : 上で得られた解を、-3 ≤ x ≤ 1.5 の範囲でプロットしてみよう
    注意
    ・微分を表すプライム記号 「'」 は単変数関数の場合のみに使用可．
    ・多変数関数の(偏)微分値を初期値として指定する場合は
    In[]:=
    D[f[x,y],x]/.x0​​Derivative[1][f][0]​​Derivative[1,0][f][0,y]
    などの表式を用いる(この3つはすべて同じ意味)．
    上の減衰振動の例で試してみると
    DSolve[{y''[t]+1/10y'[t]+y[t]0,y[0]1,y'[0]0},y[t],t]
    DSolve[{y''[t]+1/10y'[t]+y[t]0,y[0]1,(D[y[t],t]/.t->0)0},y[t],t]
    DSolve[{y''[t]+1/10y'[t]+y[t]0,y[0]1,Derivative[1][y][0]0},y[t],t]
    例題2 : x’(t) = a x(t) - b x(t)^2, x(0) = 0.1. たとえば，a=3, b=1などとして0≤t≤3でプロット．(Logistic equation)

    連立微分方程式

    xをパラメータとする関数y(x), z(x)に関する微分方程式(aは定数)
    を, 初期条件
    のもとで解こう.これらの式と初期条件はまとめて{ }の中に入れ, 解くべき関数y, zも{ }の中に入れる.
    In[]:=
    sol=DSolveValue[{
    ′
    y
    [x]y[x]-az[x],
    ′
    z
    [x]y[x]-z[x],y[0]1,z[0]4},{y[x],z[x]},x]
    In[]:=
    ParametricPlot[sol/.a2,{x,0,10}]
    ​
    Slide
    4
    of
    18

    漸近近似

    In[]:=
    AsymptoticDSolveValue[{y''[x]-y[x]0,y[0]0,y'[0]1},y[x],x0]
    AsymptoticDSolveValue[{y''[x]+y[x]0,y[0]1,y'[0]0},y[x],{x,0,10}]
    典型的には，厳密解が得られない場合や，解釈・比較のためにシンプルな表式を得たい場合に用いられる
    In[]:=
    AsymptoticDSolveValuey''[x]-y'[x]+
    10
    x
    y[x]1,y[0]1,y'[0]-34,y[x],{x,0,5}
    次数1/2のBesselの微分方程式をx=0のまわりで解く
    In[]:=
    sol=AsymptoticDSolveValue[{x^2y''[x]+xy'[x]+(x^2-1/4)y[x]0},y[x],{x,0,24}]
    In[]:=
    Plot[{(sol/.{
    
    1
    1,
    
    2
    0}),(sol/.{
    
    1
    0,
    
    2
    1})},{x,0,3π},PlotRange{-2,2}]

    クレロー(Clairaut)の微分方程式

    In[]:=
    eqn=y'[x]^2-xy'[x]+y[x]==0;
    In[]:=
    sol=DSolveValue[eqn,y[x],x]
    In[]:=
    lines=Table[cx-c^2,{c,-20,20,1/2}];​​Show[Plot[lines,{x,-10,10},PlotRange{-10,20}],Plot[x^2/4,{x,-10,10},PlotStyle{Red,Thickness[0.01]}]]
    包絡線の式を使って得ることは可能
    In[]:=
    sol/.Solve[D[sol,C[1]]==0,C[1]][[1]]

    積分微分方程式

    l=5;(*inductance*)​​r=0.5;(*resistance*)​​c=0.5;(*capacitance*)​​em[t_]=2Cos[t];(*electromotiveforce*)
    ieqn=li'[t]+ri[t]+1/cIntegrate[i[tt],{tt,0,t}]==em[t];​​ic=i[0]1;
    ​
    Slide
    5
    of
    18

    NDSolve / NDSolveValueで微分方程式を解く(DSolveでは解けない問題)

    得られた数値解を微分したり積分したりすることもできる
    例題3: 微分方程式 y'' (x) = -y (x) sin (y (x)) を, 初期条件 y (0) = 1, y' (0) = 0 の下で0≤x≤30の範囲で解いて、プロットしてみよう. これは解析的には解けない.

    連立微分方程式(NDSolve版)

    Lorentz model

    例題: 空気抵抗を受けて一様重力場中を運動する質点 (速度に比例する粘性抵抗を仮定)

    オプション (WhenEvent) を使って，地面に到達したら積分ストップ
    WhenEvent例
    (参考) 同じくWhenEventを使って，地面で跳ね返らせてみる

    2 次元正方格子中の荷電粒子の運動(WhenEvent利用．ただし隣接格子による場は無視)

    ​
    Slide
    6
    of
    18

    Schrödinger方程式(調和振動子の固有関数)

    DEigensystemは線形微分作用素の固有値，固有関数を得る

    連立微分方程式を行列表示のまま解けることの確認

    2 準位系 (行列表示)

    行列h, ベクトルψで表記

    規格化確認
    ​
    Slide
    7
    of
    18

    Shooting法

    ある種の微分方程式は同じ条件下でも複数の解をもち得る．
    次の解も上と同じ境界条件を満たす
    このような解をみつけるにはShooting法が有効．
    境界条件x(0)=x(10)=0の下で，x’(0)=1.5, 1.75, 2をShooting法のための初期値(seeds)として指定し，同じ境界条件を満たす3つの独立な解を得る
    NDEigensystemは (いまのところ) 線形のシステムにしか適用できないため，この例には使えない．

    硬い微分方程式(Stiff ODEs)

    (刻み幅をきわめて小さくしない限り，解が数値的に不安定なもの)
    硬い微分方程式として有名な例 : 化学反応を記述するRobertsonの式
    Implicit Runge-Kutta法を使うと解けるが，遅い
    偏微分方程式においても同様

    (NDSolveオプションについてのメモ) MaxSteps と MaxStepSize

    現在のMathematicaでは，多くの場合どちらもあまり大きな改善効果はなく個別の問題に大きく依存するが，問題によっては試してみる価値はある．
    MaxSteps: 最大ステップ数(刻みの数).
    Mathematicaは最初の数ステップで自動的に何ステップ必要かを推定，必要メモリ容量や計算時間を概算して許容量を超えると判断した場合は計算をあきらめる(Abortする)．ステップの幅は適応的(adaptive)に変化させるため，(ステップ幅)=(積分範囲)/ステップ数となるわけではない．
    MaxStepSize: 最大ステップ幅(刻み幅)．
    ステップの幅の最大値を与えることで，計算時間を犠牲にしてでも細かい刻みで積分していく．MaxStepSize->Infinityで制限なく細かくする．
    時刻πで撃力 (のような力)が加わった調和振動子
    ​
    Slide
    8
    of
    18

    偏微分方程式

    1次元熱伝導方程式(両端が一定温度)

    非線形偏微分方程式(KdV方程式)

    例題3: xmin = -8 から xmax = 8の範囲で，KdV方程式をt=-1からt=1において解いてみよう．初期条件はu(x,0)=2 sech(x)^2, (周期的)境界条件u(xmin,t)=x(xmax,t)
    (参考)よく見るふたつのソリトンの重ね合わせ

    周期的境界条件(波動方程式)

    境界条件の矛盾に注意

    熱伝導方程式
    おまけ: 解析解

    一見わかりにくい矛盾

    波動方程式
    境界条件・初期条件
    ​
    Slide
    9
    of
    18

    矩形領域内での波動方程式 (固定端)

    パラパラマンガ用画像を80枚作成 (要約80秒)

    アニメ化 (右のボタンでスピード調整)

    開放端

    (参考)計算の進捗モニター

    任意領域上の波動方程式

    次のセクションで！
    ​
    Slide
    10
    of
    18

    有限要素法による(非線形)偏微分方程式の求解

    有限要素法で解くことのできる偏微分方程式の標準形(“Coefficient form”)

    uはスカラー関数でもベクトル関数でもOK. (uがベクトルの場合は係数は行列)
    Navier-Stokes(非圧縮性)
    を得る

    境界条件

    ディリクレ条件
    DirichletCondition[ ]
    一般化ノイマン条件(Robin boundary condition)
    NeumannValue[ ]

    有限要素法で解くことのできる非線形偏微分方程式の係数の依存性

    ​
    Slide
    11
    of
    18
    有限要素法パッケージをロード

    Dirichlet Condition

    境界で解となる関数のとる値を指定
    predicateはeqnで指定する値が有効となる(境界上の)領域を示す条件(e.g., x≥0, x≤0&&Abs[y]<1, など)

    単位円上のポアソン方程式

    面積最小化問題

    DirichletCondition[eqn, pred] のpredをTrueとすると，境界のすべての点に適用される

    任意領域上の波動方程式

    詳細はyoutubeで解説!
    https://youtu.be/4ywRIwkAV0E
    ​
    Slide
    12
    of
    18

    Neumann Boundary Conditions

    ノイマン境界条件は領域境界に垂直な方向へ流れる流束の量を指定する
    ざっくり言って境界における解の微分値
    "標準形"の空間に関する部分(時間微分以外)を考えてみよう
    これにテスト関数ϕをかけ，領域Ωにおいて積分する
    部分積分して弱形式を得る
    と与えることの由来である．NeumannValueはDirichletConditionと違い，微分方程式の一部として指定することに注意

    NeumannValue例

    1D熱伝導方程式(断熱条件)

    1D波動方程式(吸収端)

    先ほどの2Dの例
    ​
    Slide
    13
    of
    18

    Inactive[ ] の使用

    この作用素をそのままNDSolveに与えると，標準形として解釈されない．
    ∇などの微分作用素にInactive[ ]を使って，NDSolveが式を解釈する前に式の形が変形されないようにする
    ​
    Slide
    14
    of
    18

    周期的境界条件

    1D波動方程式(再び！)

    Schrödinger方程式(2D正方形領域中の調和振動子ポテンシャル + 周期的境界条件)

    ドーナツ形の一部の上のポアソン方程式(バウムクーヘンのかけら形) + 周期的境界条件

    FindGeometricTransform[pts1, pts2]
    点列pts2をレファレンス点列pts1に移す幾何的変換を求める
    上図の青い斜めの線分を，水平の赤い線分へ移す変換
    PeriodicBoundaryConditionを設定
    ポアソン方程式
    バウムクーヘン片をつなげてリング状にして図示
    ​
    Slide
    15
    of
    18

    非線形偏微分方程式

    円柱のまわりの流れ

    Navier-Stokes方程式(非圧縮性流体)
    ​
    Slide
    16
    of
    18

    熱による密度変化(Boussinesq近似)

    PDE

    境界条件・初期条件

    時刻t=0で底面を高温，上表面を低温に設定

    NDSolve

    H. Elman, D. Silverser, and A. Wathen, “Finite Elements and Fast Iterative Sovlers”, Chapter 11
    ​
    Slide
    17
    of
    18
    ◼
  • DSolveの基本: tutorial/Calculus#9252
    • ◼
  • DSolve monograph: tutorial/DSolveOverview
    • ◼
  • PDE models monograph (PMLについての記述も): PDEModels/tutorial/Acoustics/AcousticsTimeDomain
    • ◼
  • PDEModels (自然科学で頻出するPDEモデルのコンポーネント，物質データなども利用可): PDEModels/tutorial/PDEModelsOverview
    • ​
    琵琶湖形状内で波動方程式を解いてアニメにする例の詳細はyoutubeで解説しています
    https://youtu.be/4ywRIwkAV0E
    ​
    Slide
    18
    of
    18