読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

一次元の拡散方程式をノイマン境界条件の下で解く(差分法、陰解法)

R 微分方程式 数値計算

陰解法

以下の偏微分方程式が拡散方程式である。

f:id:abrahamcow:20160816085518p:plain

ここで x は空間方向の位置を表し、t は時間を表す変数である。u(x, t) が求めたい未知関数で、初期条件 u_0(x) は与えられた関数である。D は正の定数である。

条件  u_x(t,a) = 0u_x(t,b) = 0 のように端点における微分の値を固定する条件を、ノイマン境界条件と言う。いま、この条件は境界において物の出入りがないことを意味している。

以下では拡散方程式の数値解法について述べる。

時間方向 t を刻み幅 \Delta t で分割し。上付きの添え字を用いて  t^m = m\Delta t とする。また、空間方向 x を刻み幅 ∆x = (b − a)/N で分割し、x_n = a+n \Delta x とする。

u_t = D u_{xx} を差分で近似して、

\displaystyle \frac{U_{n}^{m+1} - U_{n}^{m}}{\Delta t} = D \frac{U_{n+1}^{m}-2U_{n}^{m}+U_{n-1}^{m}}{\Delta x^2}

とする。

\displaystyle r = D \frac{\Delta t}{\Delta x^2}

とおき、拡散方程式を近似する差分方程式として、

\displaystyle (1+2r)U^{m+1}_{n} - r(U^{m+1}_{n+1}+U^{m+1}_{n-1}) = U^{m}_{n} \tag{1}

が導かれる。

問題は境界条件の扱いである。

 u_x(t,a) = 0u_x(t,b) = 0 をそれぞれ、

\displaystyle \frac{U^{m+1}_{1}-U^{m+1}_{0}}{\Delta x} = 0

\displaystyle\frac{U^{m+1}_{N}-U^{m+1}_{N-1}}{\Delta x} = 0

と近似して、

\displaystyle U^{m+1}_{0} = U^{m+1}_{1}

\displaystyle U^{m+1}_{N} = U^{m+1}_{N-1}

という差分方程式を採用することにする。

(1) 式で n = 1 とした式に、U^{m+1}_{0} = U^{m+1}_{1} を代入して、

\displaystyle (1+r)U^{m+1}_{1} - rU^{m+1}_{2} = U^{m}_{1}

同様に (1) 式で n = N − 1 とした式に、U^{m+1}_{N} = U^{m+1}_{N-1} を代入して、

\displaystyle (1+r)U^{m+1}_{N-1} - rU^{m+1}_{N-2} = U^{m}_{N-1}

こうして、N − 1 個の未知数 \{U^{m+1}_{n}\}_{i=1,2,\ldots,N−1} に関する N − 1 個の連立 1 次方程式を得る。

ベクトルと行列を用いて表わすと

\displaystyle AU = b

となる。

ただし

f:id:abrahamcow:20160816085659p:plain

f:id:abrahamcow:20160816162552p:plain

f:id:abrahamcow:20160816162543p:plain

陰解法のプログラム


R を用いる。

f:id:abrahamcow:20141124141446p:plain

溶質が均一に分布している材料 A から、材料 B に溶質が供給されていく状況を考える。

初期条件を以下のように与える。

初期条件:
\displaystyle \left\{ \begin{array}{ll} u = 1 & (x \le 10) \\ u = 0 & (x > 10) \end{array} \right.

library(tidyr)
library(cowplot)
diffeq <- function(dx, dt, tn, d, u0){ 
  n <- length(u0)
  m <- length(tn)
  r <- (d*dt)/(dx^2)

  A<-diag(1+2*r,n,n)
  for(i in 1:(n-1)){
    A[i,i+1] <- -r
    A[i+1,i] <- -r
  }
  A[1,1] <- 1+r
  A[nrow(A),ncol(A)] <- 1+r
  A_inv <- solve(A)

  b <- matrix(,m,n)
  b0 <- u0
  for(i in 1:m){
    b1 <- A_inv%*%b0
    b0 <- b1
    b[i,] <- b0
  }
  b
}
tn=1:100
u0=c(rep(1,10),rep(0,10))
out <-diffeq(dx=1,dt=1,tn=tn,d=2,u0=u0)
outdf <-as.data.frame(out)
colnames(outdf) <- seq(length.out = length(u0),by=1)
outdf2 <-data.frame(time=tn,gather(outdf, x,density))
outdf2$x <- as.numeric(outdf2$x)
ggplot(outdf2,aes(x=x,y=density,colour=time,group=time))+
  geom_line()

f:id:abrahamcow:20160816085912p:plain

付録

TeX メモ

\begin{cases}
u_t = D u_{xx} & (t>0,a<x<b)  \\
u(0,x) = u_0(x) & (a \le x \le b)\\ 
u_x(t,a) = 0 & (t>0)\\
u_x(t,b) = 0 & (t>0)\\
\end{cases}
A=\begin{pmatrix}
1+r & -r & 0 & 0 & \cdots & 0 & 0 & 0\\
-r & 1+2r & -r & 0 & \cdots &0 & 0 & 0\\
0 & -r & 1+2r & -r & \cdots &0 & 0 & 0\\
\vdots &\vdots &\vdots &\vdots & \ddots & \vdots &\vdots &\vdots\\
0 & 0 & 0 & 0 & \cdots & 1+2r & -r & 0\\
0 & 0 & 0 & 0 & \cdots & -r & 1+2r & -r\\
0 & 0 & 0 & 0 & \cdots & 0 & -r & 1+r
\end{pmatrix}
U=\begin{pmatrix}
U_{1}^{m+1} \\
U_{2}^{m+1}\\
U_{3}^{m+1}\\
\vdots \\
U_{N-3}^{m+1} \\
U_{N-2}^{m+1}\\
U_{N-1}^{m+1}\\
\end{pmatrix}
b=\begin{pmatrix}
U_{1}^{m} \\
U_{2}^{m}\\
U_{3}^{m}\\
\vdots \\
U_{N-3}^{m} \\
U_{N-2}^{m}\\
U_{N-1}^{m}\\
\end{pmatrix}

はじめまして★

ポリシー

f:id:abrahamcow:20160816024559j:plain

こんにちは~ 牛です

最近ブログのアクセス数が増えてきて、リアルのお知り合い以外の方も見てくれてるみたいでびっくりしております><

本当にありがとうございます(;_;)

さてさて
そういえば自己紹介してなかったなーなんて思いました!!!

最近まったく出会いがないなぁと思って、思い切ってはてなブログを始めました。

よく癒し系と言われます。

めっちゃ甘えたがりなので、ちゃんと受け止めてくれる人が好きです。

一緒にいるときは、ずーっとくっついていたいです(*´艸`*)

ここで運命の人に出会えるかもなんて思ってる私は甘いのかな。。

でも早く好きな人をひとりだけ見つけて退会したいです(;_;)/~~~

まずはメールから仲良くなれたら嬉しいです。よろしくお願いします(o(´∀`)o)ワクワク

女の人からの連絡も歓迎ですよ〜笑

血液型:B
星座:さそり座
興味あること:恋人、結婚相手、ドライブ、お茶したい、グルメ、飲み友、ゲーム、カラオケ、H目的お断り