image-of-web Photo by Annie Spratt on Unsplash

境界値問題 - 偏微分方程式で記述されるもう一つの典型的な問題

前回，偏微分方程式の数値解法である method of line という手法で扱った問題は時間発展問題だった．次回以降も再び時間発展問題を扱う予定であるが，今回，偏微分方程式で記述されるもう一つの典型的な問題である境界値問題を扱っておこう．

境界値問題とは，多くは偏微分方程式で記述される問題で，「境界条件を満たすような偏微分方程式の解を探す」ものである．

初期値問題も偏微分方程式と境界条件から構成されるが，それとの違いは以下のように考えられる:

初期値問題は境界条件が与えられている境界から与えられていない境界の方向へ解の依存関係の方向性がある．よって，この方向に沿って計算を進める解法を採用でき，少しずつ近似解を求めることができる．
境界値問題の解にはそのような依存関係の方向性がない．そのため，本来は制約を満たす解全体すべてを一気に求める必要がある．

であるので，境界値問題の数値計算はそのままではメモリも計算量も負荷が大きい．そのため，実際には「一度に解かずに済むように工夫した」手法を適用することが多い．まあ今回はそのあたりにあまり細かく触れずに解説しよう．

具体的問題: Poisson 問題

既知の実数値関数 $f(\boldsymbol{x})$ によって未知の実数値関数 $u(\boldsymbol{x})$ が満たすべき方程式として

\[ \Delta u(\boldsymbol{x}) = f(\boldsymbol{x}) \mbox{ in } \Omega \subset \boldsymbol{R}^n \]

があたえられるとき，この方程式を n次元 Poisson (ポアソン)方程式と呼ぶ，のだな．ただし，$\Delta$ はラプラシアンで，$\boldsymbol{x} = (x_1, x_2, \cdots, x_n)$ と書いたとき

\[ \Delta := \frac{\partial^2}{\partial x_1^2} + \frac{\partial^2}{\partial x_2^2} + \cdots + \frac{\partial^2}{\partial x_n^2} \]

のことと思って良い¹．

ただ，Poisson方程式は 2階の微分方程式だから，それだけでは関数 $u$ の値は定まらない．未知の定数が独立変数ごとに２つ残るので，それらを定めるような境界条件が必要だ．

これらの境界条件も合わせて指定されたとき，解 $u(\boldsymbol{x})$ を定める問題を Poisson問題²といい，これは典型的な境界問題である．

1次元 Poisson問題

例題として，例えば次のような1次元 Poisson問題を考えてみよう(1次元なので偏微分方程式ではなくて常微分方程式になってしまうが)．

\[ \left\{ \begin{array}{rclr} \displaystyle \frac{d^2 u}{dx^2} & = & \sin(x) & \mbox{ in } [0, L], \cr u(0) & = & 3, & \cr u(L) & = & 4. & \end{array} \right. \label{eq:1dim_poisson} \]

この問題は手で解けるように設定してあるが，しばらくそのことに気づかないふりをして話を進めよう．

さて，もしも $x=L$ での $u$ の値が与えられているのではなく $x=0$ での $du/dx$ の値が与えられていれば，それは初期値問題であり， A2. 二階微分はどうよ: Newton運動方程式の内容で扱えるのだが，ここではそうではない．

というわけで，どうやって近似解を求めるかあらためて考えないといけない．それには，まずは問題を離散近似してみるのが良い．ということで，やってみよう．

まず，空間を $N$等分する．刻み幅は $\Delta x = L/N$ になるので，その刻み $x = k\Delta x$, $k = 0, 1, \cdots N$ の上での値だけを近似値として使うことにして，その近似値を $u_k \cong u(k\Delta x)$ と書くことにしよう．

次に，二階微分 $d^2/dx^2$ を二階中心差分 $\delta_k^{(2)}$ で近似することにすると，今回の問題( eq(\ref{eq:1dim_poisson}))は素直に考えると次のように離散近似されるだろう．

\[ \left\{ \begin{array}{rclr} \displaystyle \frac{u_{k+1} - 2 u_k + u_{k-1} }{\Delta x^2} & = & \sin(k\Delta x) & \mbox{for } k = 1,2,\cdots,N-1, \cr u_0 & = & 3, & \cr u_N & = & 4. & \end{array} \right. \label{eq:discr_1dim_poisson} \]

これを少し見やすく変形して，

\[ \left\{ \begin{array}{rclr} u_0 & = & 3, & \cr u_{k+1} - 2 u_k + u_{k-1} & = & \Delta x^2 \sin(k\Delta x) & \mbox{for } k = 1,2,\cdots,N-1, \cr u_N & = & 4. & \end{array} \right. \label{eq:scheme_1dim_poisson} \]

と書き換えよう．さらにこれを線形代数で扱いやすいように, 境界条件を繰り込んで下記のように直してみよう．

\[ \left\{ \begin{array}{rclr} u_{2} - 2 u_1 + 3 & = & \Delta x^2 \sin(\Delta x), & \cr u_{k+1} - 2 u_k + u_{k-1} & = & \Delta x^2 \sin(k\Delta x) & \mbox{for } k = 2,3,\cdots,N-2, \cr 4 - 2 u_{N-1} + u_{N-2} & = & \Delta x^2 \sin((N-1)\Delta x). \end{array} \right. \]

そして，未知数を左辺に，既知量を右辺に集め直すと次のようになる．

\[ \left\{ \begin{array}{rclr} u_{2} - 2 u_1 & = & -3 + \Delta x^2 \sin(\Delta x), & \cr u_{k+1} - 2 u_k + u_{k-1} & = & \Delta x^2 \sin(k\Delta x) & \mbox{for } k = 2,3,\cdots,N-2, \cr -2 u_{N-1} + u_{N-2} & = & -4 + \Delta x^2 \sin((N-1)\Delta x). \end{array} \right. \label{eq:lm_scheme_1dim_poisson} \]

これはよく見ると連立一次方程式だ．よって，これを線形代数の記号で書き直そう．まず，未知数を集めたベクトルを

\[ \boldsymbol{u} := (u_1, u_2, \cdots, u_{N-1})^{t} \]

とし，そこにかかる係数行列を下記のような $(N-1)\times(N-1)$ サイズの対称行列で定義する．

\[ A := \left( \begin{array}{cccccc} -2 & 1 & 0 & \cdots & & & 0 \cr 1 & -2 & 1 & 0 & \cdots & & 0 \cr 0 & 1 & -2 & 1 & 0 & \cdots & 0 \cr & & \ddots & \ddots & \ddots & & \cr & & & \ddots & \ddots & \ddots & \cr & & & & 1 & -2 & 1 \cr 0 & & & & 0 & 1 & -2 \cr \end{array} \right) \]

そして右辺相当を

\[ \boldsymbol{b} := \left( \begin{array}{c} -3 + \Delta x^2 \sin(\Delta x) \cr \Delta x^2 \sin(2\Delta x) \cr \vdots \cr \Delta x^2 \sin((N-2)\Delta x) \cr -4 + \Delta x^2 \sin((N-1)\Delta x) \end{array} \right) \]

とする．すると式 ($\ref{eq:lm_scheme_1dim_poisson}$) は

\[ A \boldsymbol{u} = \boldsymbol{b}
\]

と等価なのだ．だからあとはこれを解けば，未知であった $u_1, u_2, \cdots, u_{N-1}$ が求まるはずだ．

ここで，一気に連立一次方程式を解く必要があることが境界値問題の特徴である「解全体を一度に求めないといけない」ことの反映なのだ．似たような初期値問題だったら $u_0$ → $u_1$ → $u_2$ → $\cdots$ というように順番にスカラー値を求めていけばよいのだ，ということと対照的だということを頭に入れておこう．

では早速，julia で計算してみよう．まずは定数等だ．$L = 10$ としてやってみよう．

1
2
3


L = 10
N = 1000
Δx = L/N

次に，係数行列 $A$ を定義しないといけないね．なお，途中でスカスカの行列などを特別に扱うパッケージである SparseArrays を用いるが，これは標準ライブラリとして既にインストールされているはずなので特にインストール作業は不要だ．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11


using LinearAlgebra 
# 線形代数関連の機能を使いたいので．
using SparseArrays  
# スカスカの行列を特別に扱う．メモリや計算量が節約できるよ！

sd = ones(N-2)
# 斜めの 1 成分を代入するためにこのベクトルを用意して…

A = sparse( -2 * I + diagm( 1 => sd, -1 => sd) )
# 対角成分は -2I で，副対角成分は diagm で代入しているよ．
# 残りはゼロ値でメモリや計算量がもったいないので sparse 宣言をしてスカスカ扱いをするよ．

999×999 SparseMatrixCSC{Float64, Int64} with 2995 stored entries:
⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦

最後に sparse 扱いをしたので，行列の中でノンゼロな部分だけそれっぽく書いた表示が最後に見えているね．

あとは右辺をベクトル化した量 $\boldsymbol{b}$ の定義だ．

1
2
3


b = Δx^2 * [ sin(k*Δx) for k in 1:N-1 ]
b[1] += -3
b[N-1] += -4

これで材料は揃った．あとは連立一次方程式を解いて近似解を求めるだけだ． julia だとこれは簡単で，次の1行で終わる．

1

u = A \ b

999-element Vector{Float64}:
 2.990456057689051
 2.9809131153614348
 2.9713721729004883
 2.9618342299895617
 2.9523002860120537
 ⋮
 3.955103790376946
 3.964186748688378
 3.9732187037957853
 3.9821987981082403
 3.991126179220851

上のコードの中の \ はバックスラッシュで，日本語 OS の場合は ¥ キーと同等だ．画面で表示が違ったり，キーボードで該当キーが見つからなくても慌てないようにしよう．

さて，近似解をプロットしてみよう．

1
2


using Plots
plot(u)

ふむ．境界条件も無理なく満たされていそうだ，ということが見て取れる．

さて，そして今回は問題が手で解けるように設定されていたので，それを利用して誤差の様子も見てみよう．簡単なので導出は省略するが，今回の問題 $(\ref{eq:1dim_poisson})$ の厳密解は以下のようになる．

\[ u(x) = - \sin(x) + \left(\frac{1 + \sin(L)}{L}\right) x + 3 \]

というわけで誤差を測ろう．まず厳密解を定義する．

1

exact(x) = - sin(x) + ((1+sin(L))/L)*x + 3

誤差は以下の通りだ．

1

error = [ exact(k * Δx) - u[k] for k in 1:N-1 ]

999-element Vector{Float64}:
 8.786589322795635e-8
 1.7572345356597907e-7
 2.6356434812413454e-7
 3.513802466770244e-7
 4.391628212196963e-7
 ⋮
 3.324673301108305e-7
 2.651095640793244e-7
 1.981768238401571e-7
 1.3167625478871514e-7
 6.561495968782083e-8

悪くないが，思ったより微妙に誤差が大きいな，という感じだ．

さて，誤差をプロットしてみよう．

1

plot(error)

誤差の(絶対値で)大きいところをみると，解の極小点，極大点付近だということがわかる．

もとの問題との関係を考えると，素直に「二階微分値が(絶対値で)大きいところで誤差が大きい」と言えるので，この誤差は主に「二階差分計算の誤差」が原因だ，と考えてよいだろう．

Bratu 問題

より正確には Liouville-Bratu-Gelfand 問題などと呼ばれるのだが， $\Omega \subset \boldsymbol{R}^n$ での未知関数 $u(\boldsymbol{x})$ に対する次のような境界問題がある．境界値関数 $b(\boldsymbol{x})$ は既知である．

\[ \left\{ \begin{array}{rclrr} \nabla \cdot ( \kappa(\boldsymbol{x})\nabla u(\boldsymbol{x})) + C \exp( u(\boldsymbol{x}) ) & = & 0 & & \mbox{ in } \Omega , \cr u(\boldsymbol{x}) & = & b(\boldsymbol{x}) & & \mbox{ at } \partial\Omega . \end{array} \right. \]

なお，$C$ は正定数で，$\kappa(\boldsymbol{x}) > 0$ は地点 $\boldsymbol{x}$ での拡散係数のようなものだ．

Poisson問題と異なりこの問題は非線形問題なので，解くのがなかなか難しいのだ．

1次元 Bratu 問題

$\kappa(\boldsymbol{x})$ を定数 1 に固定し，考える領域を 1次元にし，境界値も両端でゼロにしたシンプルな問題を考えてみよう．次のような問題だ．

\[ \left\{ \begin{array}{rclrr} \Delta u(x)) + C \exp( u(x) ) & = & 0 & & \mbox{ in } [0, L] , \cr u(0) & = & 0, \cr u(L) & = & 0 . \end{array} \right. \]

考え方は先の Poisson 問題と同様で良いので，差分方程式を整理した格好の近似式 $(\ref{eq:lm_scheme_1dim_poisson})$ 相当の式をいきなり書いてみよう．ただし，未知量を左辺に集めるので少し見かけが変わるぞ．

\[ \left\{ \begin{array}{rclr} u_{2} - 2 u_1 + \Delta x^2 C \exp(u_1) & = & 0 , & \cr u_{k+1} - 2 u_k + u_{k-1} + \Delta x^2 C \exp(u_k) & = & 0 & \mbox{for } k = 2,3,\cdots,N-2, \cr -2 u_{N-1} + u_{N-2} + \Delta x^2 C \exp(u_{N-1}) & = & 0 . \end{array} \right. \label{eq:lm_scheme_1dim_bratu} \]

そこで，線形代数の記号を使って書くことにすると，$A$ や $\boldsymbol{u}$ は先の Poisson問題のものと同じで良いので， $\Delta x^2 C \exp(u)$ 相当であるベクトル値関数を

\[ \boldsymbol{e}(\boldsymbol{u}) := \Delta x^2 C \left( \begin{array}{c} \exp(u_1) \cr \exp(u_2) \cr \vdots \cr \exp(u_{N-2}) \cr \exp(u_{N-1}) \end{array} \right) \]

とする．すると式($\ref{eq:lm_scheme_1dim_bratu}$)は次のように書ける．

\[ A \boldsymbol{u} + \boldsymbol{e}(\boldsymbol{u}) = \boldsymbol{0} \label{eq:abstract_lm_scheme_1dim_bratu} \]

これは連立非線形方程式なので，NLsolve package の力を借りて解くことになるだろうね．

では Julia で実行してみよう．なお，定数 $C$ は今回は $C = 1.0$ として解いてみることにしよう．

まず，定数等や行列の設定は先の Poisson 問題と同様で良いが，非線形問題なのでサイズは小さくしておこう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11


C = 1.0
 
L = 1    # 領域を小さく
N = 100  # 変数の数も少なく
Δx = L/N

using LinearAlgebra 
using SparseArrays  

sd = ones(N-2)
A = sparse( -2 * I + diagm( 1 => sd, -1 => sd) )

次に関数 $\boldsymbol{e}$ を定義しよう．

1

e(u) = Δx^2 * C * exp.(u)

そして，ゼロになってほしい，と NLsolve に渡す関数だ．まあ，式($\ref{eq:abstract_lm_scheme_1dim_bratu}$) のことだが．

1

f_zero(u) = A*u + e(u)

ではいつものように NLsolve を使う準備をしよう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12


using NLsolve

function nls(func, params...; ini = [0.0])
    if typeof(ini) <: Number
        r = nlsolve((vout,vin)->vout[1]=func(vin[1],params...), [ini])
        v = r.zero[1]
    else
        r = nlsolve((vout,vin)->vout .= func(vin,params...), ini)
        v = r.zero
    end
    return v, r.f_converged
end

次に，NLsolve に渡す初期値を適当に作ろう．まあ，境界条件を満たしていて，かつ，定数ではないほうが良いだろうから，例えば

\[ u_0(x) = 0.1 \left\{ \frac{L^2}{4} - \left(x- \frac{L}{2} \right)^2 \right\} \]

などの感じで良いのではないかな．ではこれをプログラムだ．

1

u_0 = [ 0.1 * (L^2/4 - (k*Δx-L/2)^2 ) for k in 1:N-1 ]

99-element Vector{Float64}:
 0.000990000000000002
 0.001960000000000001
 0.0029100000000000016
 0.0038399999999999992
 0.004749999999999999
 ⋮
 0.004749999999999994
 0.003840000000000005
 0.0029100000000000016
 0.001960000000000001
 0.000990000000000002

では解いてみよう．

1
2


u, flag = nls(f_zero, ini = u_0)
# 成功したかどうか(最後に true と出るかどうか)を見るためこのようにしている．

([0.005443486690315003, 0.010786427547691988, 0.016028283923992016,
 0.021168524557703402, 0.026206625774707157, 0.03114207168783226, 
 0.03597435439505224, 0.040702974176173724, 0.04532743968786691,
 0.04984726815688744  …  0.049847268156888044, 0.04532743968786747, 
 0.040702974176174224, 0.03597435439505269, 0.03114207168783265, 
 0.026206625774707483, 0.021168524557703662, 0.016028283923992213, 
 0.01078642754769212, 0.005443486690315069], true)

"true" と最後に表示されているので計算に成功したようだ．では早速プロットしてみよう．

1
2


using Plots
plot(u)

ふむ，なるほど，こんな感じか．

さて，ここで少し考えたい．そもそもこの解はこの問題の唯一の解だろうか．実はそうではないのだ．実は $C = 1.0$ のときはもう一つ解があって，これよりだいぶ大きいのだ．

そこでそれも求めてみよう．それには，NLsolve に渡す初期値を例えば以下のように変えてみれば良い．ただし，こうしたやり方はいつもうまくいくとは限らない(というか，失敗することのほうが多いかな…).

1
2


u_0 = [ 20 * (L^2/4 - (k*Δx-L/2)^2 ) for k in 1:N-1 ]
# さっきの初期値の200倍だ．

早速解いてみよう．

1

u, flag = nls(f_zero, ini = u_0)

([0.10840242607728252, 0.21669340253881816, 0.3248601826742571, 
0.43288857909489387, 0.5407628050722325, 0.6484652994155491, 
0.755976533428671, 0.8632747984196193, 0.9703359721851871, 
1.077133262862372  …  1.0771332628623853, 0.9703359721851992, 
0.8632747984196301, 0.7559765334286803, 0.6484652994155571, 
0.5407628050722392, 0.4328885790948992, 0.3248601826742611, 
0.2166934025388208, 0.10840242607728384], true)

うまく計算できたことを確認して，プロットだ．

1

plot(u)

どうやらうまくいったようだな．

簡単な印象だが？

…と手始めにやってみたが，どうも簡単な印象を受ける．境界値問題は本当に簡単なのだろうか？

実はそうでもないのだ．これまでのものは「意図的に簡単な例を挙げていた」のだ．

というわけで，問題そのものを変えずに難しさを少しだけ味わってみよう．例えば，上の 1次元 Bratu 問題で定数 $C$ の値を

\[ C = 3.5137 \]

にして，同様に問題を解いてみよう(このとき解があることは実はわかっている)．

1
2
3
4
5
6


C = 3.5137

u_0 = [ 0.1 * (L^2/4 - (k*Δx-L/2)^2 ) for k in 1:N-1 ]
# u_0 については自分の好みで変えて良いぞ．

u, flag = nls(f_zero, ini = u_0)

([0.039819444888453855, 0.07927324615921884, 0.11834668946240666, 
  0.1570246187738921, 0.19529143679277242, 0.23313110677508755, 
  0.27052715595947624, 0.3074626806699197, 0.343920353270724, 
  0.37988243104433617  …  0.37988243104438174, 0.34392035327115206, 
  0.30746268066992904, 0.27052715595979726, 0.23313110677507384, 
  0.19529143679299493, 0.15702461877386875, 0.11834668946254003, 
  0.0792732461591979, 0.039819444888507534], false)

う… 計算に失敗する．

さてこうなると実感できるが，境界値問題に対して他の解法を知らないとここで手詰まりになりかねない．結構厄介なのだ．

自分で解法を工夫しても良いし，他の解法を調べてそれを使っても良いので，この $C$ のもとで 1次元 Bratu問題に対してなんとか近似解を求めてみよう → レポート問題

(おまけ) 単純な 1次元 Bratu 問題についての知見

実はこのケースについては少し調べがついていて，係数 $C$ について臨界値 $C_c := 3.513830719\ldots$ が存在して

$C < C_c$ の場合: 実数解関数 $u(x)$ が 2つ存在する．
$C = C_c$ の場合: 実数解関数 $u(x)$ が 1つ存在する．
$C > C_c$ の場合: 複素数値解関数 $u(x)$ が 2つ存在する．

ということがわかっている．しかも解の形もわかっていて，$C$ に対して

\[ \cosh(\alpha) = \frac{4}{\sqrt{2C}}\alpha \]

を満たす $\alpha$ (これが複数あったりするのだ)を用いて

\[ u(x) = 2 \log\left( \frac{ \cosh(\alpha) }{ \cosh( \alpha(1-2x) )} \right) \]

と書けるのである．だから誤差の計算も可能だね．

レポート

下記要領でレポートを出してみよう．

e-mail にて，
題名を 2021-numerical-analysis-report-10 として,
教官宛(アドレスは web の "TOP" を見よう)に，
自分の学籍番号と名前を必ず書き込んで，
内容はテキストでも良いし，pdf などの電子ファイルを添付しても良いので，

下記の内容を実行して，結果や解析，感想等をレポートとして提出しよう．

1次元 Poisson問題の右辺を $\cosh(3x)$ にして数値計算を実行してみよう．できれば誤差もプロットしたい．
ただし，$x \in [0,10]$ とし，境界条件は $u(0) = 3, u(10) = 4$ とする．
(チャレンジ問題) どのような方法でも良いので，1次元 Bratu問題の近似解を $C = 3.5137$ の条件下でなんとかして数値計算で求めよう．
(チャレンジ問題) 境界値問題に対する Shooting method を調べて，今回の問題のいずれかに適用してみよう．
(チャレンジ問題) 2次元 Poisson問題(右辺は適当に設定して良い)を用意して，数値計算してみよう．数式云々よりもプログラムが厄介だろうからそこに留意しよう．

より正確には，「空間が平坦なときは」という前提条件付きだ． ↩︎
実は，Poisson 問題には Green関数と呼ばれる特殊な関数 $G(x,s)$ を用いた基本解(境界値をゼロにしたときの厳密解相当) $u_\mbox{b}(x) = \int_{\Omega} f(s)G(x,s)ds$ を求め，あとは境界条件を満たすように一次関数相当を足すだけ，という方法がある．1次元の場合は Heaviside の階段関数を $H(t)$ として $G(x,s) = (x-s) H(x-s)$ なので，$u_\mbox{b}(x) = \int_\mbox{leftside}^x f(s)(x-s)ds$ になる． ↩︎

数値計算法基礎 (2021)

10. 境界値問題 via 偏微分方程式