So-net無料ブログ作成
数値解析 ブログトップ
前の10件 | -

超越方程式の近似解法 [数値解析]

超越方程式の近似解法

 

f(x)は何回でも微分可能な関数とする。

方程式

  

の一つの解をαとすると、

  

さて、点x₀でテーラー展開すると

  chokin-001.png

2次の項を無視し、さらに、f(x)=0とすると、

  

という近似式が得られる。

Chokin-Graph-000.png幾何的に言うと、(x₀,f(x₀))におけるy=f(x)の接線の方程式は

  

だから、接線のx切片は

  

つまり、関数f(x)x₀におけるf(x)の接線で近似し、その接線のx切片を方程式f(x)=0の近似値としている。

 

(1)式で得られたxx₁とし、

  

と、さらに

  chokin-000.png

と逐次的に計算し、方程式f(x)=0の解αを数値的に求める方法をニュートン法ニュートン・ラプソン法)という。

 

問1 ニュートン法を用いて、01の間にある

  

の(近似)解を求めよ。

【解】

  

x₀=0とすると、(2)式より

  

x₁=0.5だから、

  

x₂=0.61だから

  

よって、x=0.619が解である。

(解答終)

 

このように、ニュートン法は数回の計算で方程式f(x)=0の近似解を求めることができる。

 

しかし、一々、を計算することが面倒な場合、(2)式ではなく

  

を用いて、f(x)=0の近似解を求めることもできる。

この方法をvon Misesという。

 

問2 von Mises法を用いて、01の間にある

  

の(近似)解を求めよ。

【解】

  

x₀=0とすると、

  

したがって、(3)式より

x₁=0.5とすると、

  

以下、同様に

  chokin-002.png

よって、x=0.619

(解答終)

 

上の問からわかるように、一般に、von Mises法はニュートン法よりも収束速度は遅い。

 

下図のように、方程式f(x)=0の解αの間にあるとする。

 

Chokin-Graph-001.png

 

このときを結ぶ直線の方程式は

  

このx切片を求めると、

  chokin-003.png

とすると、

  

という漸化式が得られる。

このように、1次補間を用いて、f(x)=0の(近似)解を求めることもできる。

 

問3 (4)を用いて、01の間にある

  

の(近似)解を求めよ。

【解】

x₀=0x₁=1として表計算ソフトを用いて計算した結果を以下に示す。

 

chokin-tab-002.png

 

(解答終)

 

上の表のn=3のとき(黄色の行)、解は0x₃≒0.497の間に存在せず、そのため、n=4(水色の行)のとき、外挿によって方程式の解の近似値が計算されていることに注意。x₀x₁の選び方によってはこのようなことが起きる。

ニュートン法ほど速くないけれど、それでも、5〜6回ほどの計算で収束していることがわかる。

 

 


nice!(0)  コメント(0) 

境界値問題の前進積分による解法2 [数値解析]

境界値問題の前進積分による解法2

 

次の微分方程式があるとする。

  

この解は、

  

 

2階常微分方程式(1)の境界値問題を4次のルンゲ・クッタ法で解く方法を思いついたので、紹介する。

微分方程式

  

は、

  

とおくと、

  

となり、次の連立微分方程式に書き換えることが可能。

  

そして、(2)は、x=x₀におけるypの初期値y₀p₀が与えられれば、(2)は(4次の)ルンゲ・クッタ法を用いて前進的に解くことができる。

しかし、(1)の境界値問題では、pの初期値p₀が与えられていない。

そこで、y₀=0p₀=1と推測し、n=10h=0.1として、ルンゲ・クッタ法で解くと、次のような計算結果が得られる。

 

RunRun-tab-001.png

 

x=1におけるyの値はy(0)=1だから、として計算したyの誤差

  

次に、とし、ルンゲ・クッタ法を用いて解くと、次のようになる。

 

RunRun-tab-002.png

 

このとき、x=1におけるyの計算値は1.425591だから、誤差

  

である。

この計算結果を用いてx=0におけるpを次のように推測する。

  

になる。

ところで、この微分方程式の厳密解は

  

であり、これを用いてy'(0)を求めると、

  

だから、何と、小数点5位まで正確な値を推測できる。

先に求めたを小数点7位で四捨五入した

  

として計算した結果は次の通り。

 

RunRun-tab-003.png

 

x=1でのyの計算値は1になっており、この境界値問題を(4次の)ルンゲ・クッタ法を用いて見事解くことができた!!

 

3回も計算するのは面倒なので、推測機能を備えたプログラムを新たにC言語で作った。

 ――最近、Fortranでずっとプログラムを書いていたから、Cでプログラムを書くことに、結構、手こずった(^^ゞ――

それを用いて解いた結果は次の通り。

 

RunRun-tab-004.png

 

計算に使用したプログラムは次の通り。

 

#include <stdio.h>
#include <math.h>
#define N 10

double f(double , double, double);
double g(double , double, double);
double exact(double);

void Runge_Kutta(double *x, double *y, double *z, double h, int n) {
    double dk[2][4];
    int i;
   
    for (i=0; i < n; i++) {
        dk[0][0] = h*f(x[i],y[i],z[i]);
        dk[1][0] = h*g(x[i],y[i],z[i]);
        dk[0][1] = h*f(x[i]+h/2, y[i]+dk[0][0]/2, z[i]+dk[1][0]/2.);
        dk[1][1] = h*g(x[i]+h/2, y[i]+dk[0][0]/2, z[i]+dk[1][0]/2.);
        dk[0][2] = h*f(x[i]+h/2, y[i]+dk[0][1]/2, z[i]+dk[1][1]/2.);
        dk[1][2] = h*g(x[i]+h/2, y[i]+dk[0][1]/2, z[i]+dk[1][1]/2.);
        dk[0][3] = h*f(x[i]+h, y[i]+dk[0][2], z[i]+dk[1][2]);
        dk[1][3] = h*g(x[i]+h, y[i]+dk[0][2], z[i]+dk[1][2]);
       
        x[i+1]=x[i]+h;
        y[i+1]=y[i]+(dk[0][0]+2*dk[0][1]+2.*dk[0][2]+dk[0][3])/6.;
        z[i+1]=z[i]+(dk[1][0]+2*dk[1][1]+2.*dk[1][2]+dk[1][3])/6.;
    }
}

main() {
    double x[N+1],y[N+1],p[N+1];
    double h=0.1;
    double p0=1., p1=2., w;
    double eps0, eps1, eps=0.000001;
    double yn;
    int i, cnt, n = N;


    x[0]=0.; y[0]=0.; p[0]=p0; yn=1.0; // p0はx=0におけるp=dy/dxの推測値、yn=y(1)はx=1におけるyの値
    Runge_Kutta(x, y, p, h, n);  // y0=1、p0=1としてルンゲ・クッタ法で計算
    eps0=y[n]-yn; // x=1におけるyの誤差を計算
    p[0]=p1; //p₀の推測値を更新
    Runge_Kutta(x, y, p, h, n);  // y0=1、p0=2としてルンゲ・クッタ法で計算
    eps1=y[n]-yn; // x=1におけるyの誤差を計算

    w=(eps1*p0-eps0*p1)/(eps1-eps0); // 補間法を用いてp₀を推測
    p0=p1; p1=w;
   
    cnt=1;
    while(cnt <= 20) {
        eps0=eps1;
        p[0]=p1;
        Runge_Kutta(x, y, p, h, n); // 推測値を用いて再計算
        eps1=y[n]-yn;  // 誤差を計算
        if (fabs(eps1)<eps) break; // 収束条件を満たしていたら計算終了
        w=(eps1*p0-eps0*p1)/(eps1-eps0); // 補間を用いてp₀を推測
        p0=p1; p1=w;
        cnt++;
    }

    printf("    x       y      厳密解\n");
    for(i=0;i<=n;i++) {
        printf("%f %f %f\n",x[i],y[i],exact(x[i]));
    }
}

double f(double x,double y, double z) {
    return z;
}

double g(double x, double y, double z) {
    return -2.*x*z - x*x*y;
}

double exact(double x) {
    double e;
    e=exp(1);
    return pow(e,3./2.)/(e*e-1)*(exp(x)-exp(-x))*exp(-x*x/2);
}

 

そして、「どうだ」と、胸を張るネムネコであった。

 

 

画像元:YouTubeの上の動画

 

さらに、自画自賛ソングを!!

 

 

なお、ループ(while)が存在するのは、非同次線形微分方程式だけではなく、非線形の微分方程式にも対応できるようにしたため。

線形微分方程式ならば、whileループは不要!!

 

 


nice!(0)  コメント(0) 

[2階常微分方程式] [数値解析]

2階常微分方程式]

 

 

ですか。2階常微分方程式だと境界条件(2)(初期条件)は、y'(0)0みたいになるのが普通なんですが一般論という事で(^^)

 

 まず(1)を普通に差分します。x00x1Δxx22Δx,・・・,xjjΔx,・・・,xnnΔx1として、

 

 じつは線形微分方程式と線形差分方程式は、ほとんどパラレルな形で解けます(ネコ先生はご存知でしょうが(^^))。

(3)の特性方程式は、

なので、

 

 

 よって(3)の一般解は、(5)μμ1μ2として、

 ここにC1C2は境界条件から決まる任意定数。境界条件y(0)y01より(6)から、

 境界条件y'(1)y'(nΔx)0より直接、

 (8)(6)を代入すれば、

 (7)(9)より、

 

 (6)(10)(11)を使えば、図-1が得られます。Δx0.1n10です。やっぱり、y1階微分に最後まで中

心差分を使うのが精度の決め手みたいですね(^^)

 


(執筆:ddt³さん)


nice!(2)  コメント(0) 

境界値問題の前進積分による解法1 [数値解析]

境界値問題の前進積分による解法

 

微分方程式の境界値問題を数値的に解くには、差分法などを使い離散化した連立方程式を解く必要があった。しかし、2階の境界値問題は、初期値問題のように、前進積分法によって解ける場合がある。その方法を紹介する。

 

問題 非線形境界値問題

  

h=0.2として、前進積分法で解け。

【解】

  

とすると、

  

h=0.2だから、

  

という漸化式が得られる。

問題の条件よりy(0)=0だからy₀=0

ここで、y₁=0.3と仮定すると、(2)式よりy₂

  

以下同様に、

  

が得られる。

1回目のy₁の推測値を後で使うので、と書くことにする。

y(1)=1だから、

  

y₁=0.3だと、y₅y(1)=1を超過するので、y₁0.3より少し減らし、y₁=0.25とし再計算する。

  

したがって、

  

ここで、補間を使って、y₁を次のように推測する。

  

この値を用いて、また、再計算する

  

(解答終)

 

以下に、表計算ソフトを用いて計算した結果を示す。

 

 

 

あくまで、このように計算できる場合もあるという話!!

 

ということで、これまでのように、連立方程式を解く方法でこの問題を解くことにする。

 

微分方程式を差分で近似すると、

  

いつも同じ方法だと芸がないということで、今回は、この連立2次方程式をニュートン法で解いた。

やっぱ、ニュートン法は、収束が早いね。

 

以下に、そのプログラムを示す。

 

parameter(n=10)
real y(0:n)
real a(n),b(n),c(n),d(n)
real w(n)

eps=1.e-6

! 初期化
y=0

! 境界条件
y(0)=0.; y(n)=1.0

h=1./n

! ニュートン法
do k=1,10
    do i=1,n-1
        a(i)=1.
        b(i)=2.*(h*h*y(i)-1)
        c(i)=1.
        d(i)=y(i-1)+(h*h*y(i)-2.)*y(i)+y(i+1)
    end do
   
    call tdma(a,b,c,d,n-1)

    e_max = 0.
    do i=1,n-1
        y(i)=y(i)-d(i)
        e_max=amax1(e_max,abs(d(i)))
    end do
!    write(*,*) k, e_max
    if (e_max.lt.eps) exit
end do

! 結果の出力
write(*,*) ' k      x       y'
do i=0,n
    write(*,100) i, i*h, y(i)
end do

100 format(i3,1x,f8.5,1x,f8.5)

end

!  TDMA
subroutine tdma(a,b,c,d,n)
real a(n), b(n), c(n),d(n)
do i=1,n-1
    ratio=a(i+1)/b(i)
    b(i+1)=b(i+1)-ratio*c(i)
    d(i+1)=d(i+1)-ratio*d(i)
end do
d(n)=d(n)/b(n)
do i=n-1,1,-1
    d(i)=(d(i)-c(i)*d(i+1))/b(i)
end do

end

 

計算結果は、次の通り。

 

Hisen-Result.png

 

ほとんど直線的に変化しているということもあるのだけれど、問題の解答のような粗い計算でも、小数点3位くらいまではこの数値解と一致しており、結構、いい感じで解けていることがわかる。

 

問題の微分方程式

  

は、一見、解くのは簡単そうに見えるけれど・・・。

解けると思うヒトは、是非、チャレンジして欲しい。

そして、もし、運良く解くことができたら――正確には、三角関数や指数、対数関数などの初等関数を用いて、この解を表わせたら――、教えて欲しいにゃ。

 

 


nice!(0)  コメント(0) 

差分法による2階常微分方程式の境界値問題の数値解法3 [数値解析]

差分法による2階常微分方程式の境界値問題の数値解法3

 

前回、次の微分方程式

  

x=1における第2種の境界条件(ノイマンの境界条件)

  

を後退差分

  

と近似し計算したところ、うまくいかなかった。(下図参照)

 

kyoukaichi-graph-002.png

 

これは、

(1)式を差分方程式にするときに用いた

  

であるのに対し、x=1における境界条件の近似で用いた

  

の誤差がO(h²)と、(2)、(3)式に対して誤差が大きく、これが数値解の精度を悪化させたためである。

したがって、この問題を解決するためには、x=1における境界条件の近似を(4)より高精度にする必要がある。

そこで、

  

[0,1]の外に仮想の点を設け、さらに、微分方程式

  

x=1でも成立すると仮定する。

すると、x=1における差分方程式は

   

になる。

x=1における境界条件に中心差分を用いると

  

これを(5)に代入すると、

  

これと、k=1n−1に関する差分方程式

  

を合わせると、未知数の個数はn個、対して、方程式の数もnとなるので、次の連立方程式を解くことによって、を全て求めることができる。

  

この方針に従ってプログラムを書き換え、n=10h=0.1の場合について計算した結果は次の通り。

前回と異なり、良好な計算結果が得られていることがわかるだろう。

 

Kyoukaichi-Graph-001.png

 

なお、本計算に使用したプログラムは、(1)より一般の

  

を解くものであり、以下に、計算結果とプログラムを示す。

 

Keisan-Kekka-001.png 

 

! 差分法を用いて y''+f(x)y'+g(x)y=h(x) の境界値問題(ノイマン条件)を解くプログラム
! ただし、固有値問題は解けない!!
parameter (ns=100)
real x(0:ns), y(0:ns)

x = 0.; y =0; ! 変数の初期化

n=10 ! nは分割数

a=0.; b=1. ! 境界のx座標
x(0)=a; x(n)=b
y(0)=1. ! 境界条件

call Solver(x,y,n)

! 計算結果の出力
write(*,*) ' i      x         y      Exact'
do i=0, n, 1
    write(*,100) i, x(i), y(i), exact(x(i))
end do

100 format(i3,1x,f9.5,1x,f9.5,1x,f9.5)
end

function exact(x)
e=exp(1.)
exact=2/(2-e)*exp(-x)- e/(2-e)*exp(-2.*x)
end

function f(x)
f=3.
end

function g(x)
g=2.
end

function h(x)
h=0.
end

! ここより下はいじると危険
! ブラックボックスとして使うべし


subroutine Solver(x,y,n)
real a(n),b(n),c(n),d(n)
real x(0:n), y(0:n)

n1=n-1
dx=(x(n)-x(0))/n
dx2=dx*dx

do i=1,n1
    x(i)=x(0)+i*dx
end do

! 差分法によって得られる連立方程式の係数の計算
do i=1,n1
    a(i)=1-dx/2.*f(x(i))
    b(i)=-(2-dx2*g(x(i)))
    c(i)=1+dx/2.*f(x(i))
    d(i)=dx2*h(x(i))
end do

! 境界条件
    d(1)=d(1)-a(1)*y(0) ! x=0 ディリクレ条件
!    x=1 dy/dx=0 ノイマン条件
    a(n)=2.
    b(n)=(dx2*g(x(n))-2.)
    c(n)=0.
    d(n)=dx2*h(x(n))

call tdma(a,b,c,d,n) ! TDMAで連立方程式を解く

do i=1,n
    y(i)=d(i) ! 計算結果をセット
end do

end

!  TDMA
subroutine tdma(a,b,c,d,n)
real a(n), b(n), c(n),d(n)
do i=1,n-1
    ratio=a(i+1)/b(i)
    b(i+1)=b(i+1)-ratio*c(i)
    d(i+1)=d(i+1)-ratio*d(i)
end do
d(n)=d(n)/b(n)
do i=n-1,1,-1
    d(i)=(d(i)-c(i)*d(i+1))/b(i)
end do

end

 

x=0.3近傍、y=0付近での数値解と厳密解との差が多少気になるので、n=100h=0.01で計算した結果を以下に示す。

 

Keisan-Kekka-002.png

 

h=0.1のとき、小数点2桁目に、h=0.01の場合は、小数点4桁目で、数値解と厳密解が異なっているが、数値計算の精度がこの程度だから、これはしょうがない。

 

 


nice!(1)  コメント(0) 

差分法による2階常微分方程式の境界値問題の数値解法2 [数値解析]

差分法による2階常微分方程式の境界値問題の数値解法2

 

 

まず、いきなり、簡単な微分方程式の問題!!

 

問題 次の微分方程式を解け。

【解】

微分方程式

  

の特性方程式は

  

したがって、①の基本解はなので、①の一般解は

  

 

(1) 境界条件がy(0)=1,y(1)=0なので、②より

  

これを解くと

  

したがって、

  

 

(2) ②より

  

境界条件はy(0)=1y'(1)=0だから

  

これを解くと

  

よって、

  

(解答終)

 

こんな問題を解きたいわけではなく、問題の(2)を差分法を使って解くことが今回のテーマ。

いきなり、(2)を解くのは大変なので、問題の(1)について考えることにする。

 

閉区間[0,1]n等分し、

  

さらに、

  

と書くことにする。

 

差分法を用いると、

  

と近似することができるので、微分方程式

  

の近似式は

   

になる。

未知数は、であり、連立方程式(1)の式の本数はn−1だから、連立方程式(1)を解くことによって、微分方程式を数値的に解くことができる。

 

プログラムを新たに作ってもいいけれど、過去に作ったより一般的な

  

を解くプログラムを再利用し、[0,1]n=10h=0.1と10分割し計算した結果は次の通り。

n=10h=0.1と粗い計算でも、微分方程式の解を正確に再現していることがわかるだろう。

 

 

計算に使用したプログラム

 

! 差分法を用いて y''+f(x)y'+g(x)y=h(x) の境界値問題(ディリクレ条件)を解くプログラム
! ただし、固有値問題は解けない!!
parameter (ns=100)
real x(0:ns), y(0:ns)

x = 0.; y =0;

n=10

a=0.; b=1. ! 境界のx座標
x(0)=a; x(n)=b
y(0)=1.; y(n)=0. ! 境界条件

call Solver(x,y,n)

write(6,*) ' i      x         y'
do i=0, n
    write(6,100) i, x(i), y(i), exact(x(i))
end do

100 format(i3,1x,f9.5,1x,f9.5,1x,f9.5)
end

function exact(x)
e=exp(1.)
exact=1/(1-e)*exp(-x)- e/(1-e)*exp(-2.*x)
end

function f(x)
f=3.
end

function g(x)
g=2.
end

function h(x)
h=0.
end

! ここより下はいじると危険
! ブラックボックスとして使うべし


subroutine Solver(x,y,n)
real a(n),b(n),c(n),d(n)
real x(0:n), y(0:n)

n1=n-1
dx=(x(n)-x(0))/n

do i=1,n1
    x(i)=x(0)+i*dx
end do

! 差分法によって得られる連立方程式の係数の計算
do i=1,n1
    a(i)=1-dx/2.*f(x(i))
    b(i)=-(2-dx*dx*g(x(i)))
    c(i)=1+dx/2.*f(x(i))
    d(i)=dx*dx*h(x(i))
end do
    d(1)=d(1)-a(1)*y(0) ! 境界条件
    d(n1)=d(n1)-c(n1)*y(n) ! 境界条件

call tdma(a,b,c,d,n1) ! TDMAで連立方程式を解く

do i=1,n1
    y(i)=d(i) ! 計算結果をセット
end do

end

!  TDMA
subroutine tdma(a,b,c,d,n)
real a(n), b(n), c(n),d(n)
do i=1,n-1
    ratio=a(i+1)/b(i)
    b(i+1)=b(i+1)-ratio*c(i)
    d(i+1)=d(i+1)-ratio*d(i)
end do
d(n)=d(n)/b(n)
do i=n-1,1,-1
    d(i)=(d(i)-c(i)*d(i+1))/b(i)
end do

end

 

このプログラムを利用すれば、問題の(2)の微分方程式の境界値問題も簡単に解けそうに思うだろう。

しかし、そうは問屋が卸してくれない。

(1)と(2)では、x=1における境界条件がy(1)=0y'(1)=0と境界条件の種類が違い、(2)ではの値が与えられていないからだ。

したがって、y'(1)=0からを求める方法をあらたに考えないといけない。

 

この最も簡単な解決方法は、x=1における境界条件y'(1)=0を、後退差分を用いて

  

と書き換え、(1)式に(2)式を新たに加えて解くというもの。

これで、未知数の数がn、連立方程式の数がnと一致し、解くことができるはずである。

 

この考え方に基づいてプログラムを書き換えて、n=10h=0.1の条件で解いたものが次の通り。

 

 

解くことはできたが、厳密解と数値解が一致しないので、n=100h=0.01の条件で再計算したものは次の通り。

 

 

これくらい計算格子を細かくとれば、それなりに良好な計算結果が得られるが、このような小手先の変更では、実用に足りないことがわかるだろう。

 

ちなみに、

問題の(1)のような境界条件をディリクレ条件、(2)のx=1のように1階の微分で境界条件が与えられるものをノイマン条件というが、ノイマン条件をプログラムにうまく入れるのは、結構、厄介なのだ。

そもそも、(1)式と(2)式は、必ずしも、整合していないし・・・。

 

参考までに、問題(2)のx=1における境界条件を

  

とディリクレ条件に直し、n=10h=0.1で計算したものを以下に示す。

 

 


nice!(0)  コメント(0) 

定常1次元熱伝導方程式の解法2 [数値解析]

定常1次元熱伝導方程式の解法2

 

次のように、熱伝導率が不連続に変化する場合の定常1次元の熱伝導方程式について考える。

  

x=l₁で、Tが連続であれば、x=l₁における温度T₁とすると、熱流束に関して

  

が成立し、

  

となる。

特に、l₁=l₂=l/2の場合、

  

ここで、

  

とおくと、形式的に

  

と書くことができる。

 

(3)は調和平均と呼ばれるもので、一般に

  

が成立する。

だから、

このように熱伝導率が不連続的に変化する場合、熱流束(温度勾配)の計算にあたって、熱伝導率の算術平均を使うことはできない。

 

ネムネコは前回、

  

の熱伝導率の計算において

   

と算術平均を用いたが、「この熱伝導率の計算には、算術平均ではなく調和平均を使うべきなのではないか」という話があって、この計算で算術平均を使う数学的な根拠は、結構、脆弱(^^

 

ではあるが、熱伝導率の計算に算術平均を使った場合は以下の通り。

sanjutsu.png

調和平均を使った場合は、

Chouwa.png

と算術平均を用い場合のほうが、調和平均を用いたものよりも、熱伝導率が温度に対して直線的に変化する場合、良好な計算結果を得ることができる。

熱伝導率などの物性値は、相変化や化学変化などを伴わない場合、温度変化の幅が狭いならば、ほとんど直線的に変化するので、算術平均を使ったほうがより現実に近い。

 

これで、算術平均を安心して使うことができる。

 


nice!(0)  コメント(0) 

定常1次元熱伝導方程式の数値計算 [数値解析]

定常1次元熱伝導方程式の数値計算

 

熱伝導率λが温度Tなどによって変化する場合、非定常の1次元熱伝導方程式は次のようになる。

  netsu-tei-000.png

ここで、ρは密度、cは比熱、は単位時間単位体積あたりに発生(消滅)する熱量である。

したがって、熱の発生がない場合、1次元の定常熱伝導方程式は

  

となり、境界条件が与えられれば、この解を求めることができる。

たとえば、熱伝導率λが一定で、x=0における温度がT₀x=lにおける温度がT₁の場合、

  

さらに、単位時間単位面積当たりに通過する熱量(熱流束)

  

である。

 

熱伝導率が一定の場合、温度Tは直線的に変化するので、わざわざ数値計算をする必要はない。

そこで、熱伝導率λが温度によって変化する場合の定常1次元熱伝導方程式の差分法を用いた数値解法について考えてみる。

 

(2)は

  

となるので、

  

 

の場合、差分法を用いて

  

と近似することが可能。

しかし、プログラムがすこし複雑になるので、今回、この計算法は採用しないことにする。

 

そこで、今回は、

  

という差分を用い、(2)式を

  

ここで、は、それぞれ、の中点における熱伝導率。

 

それで、

   

そして、熱伝導率λ

  

として、これを解くプログラムを作ってみた。

ちなみに、この微分方程式の解は

  

である。

 

計算領域0≦x≦1n=10分割した計算結果は、次の通り。

青い直線は熱伝導率が一定の場合。

 

 

計算に使用したプログラムは次の通り。

 

parameter (n=10)
real t(0:n), gam(0:n)
real a(n),b(n),c(n),d(n)

eps=1.e-6

dx =1./n

! 変数の初期化
t=0.; gam=1.
a=0.; b=0.; c=0.; d=0.

! 境界条件
t(n)=1.

do k=1, 10
! 熱伝導率Γを計算
    do i=1, n
        gam(i)=1+t(i)
    end do
! 連立方程式の係数を計算
    do i=1, n-1
        dx2=dx*dx
        aw=0.5*(gam(i-1)+gam(i))/dx2
        ae=0.5*(gam(i)+gam(i+1))/dx2
        ap=aw+ae
        a(i)=-aw
        b(i)=ap
        c(i)=-ae
        d(i)=0.
    end do
    a(1)=0.; d(1)=0.5*(gam(0)+gam(1))/dx2*t(0)
    c(n-1)=0.; d(n-1)=0.5*(gam(n-1)+gam(n))/dx2*t(n)

! TDMAで連立方程式を解く
    call tdma(a,b,c,d,n-1)

! 計算結果の更新と収束判定
    err=0.0
    do i=1,n-1
        err=amax1(err,abs(1-t(i)/d(i)))
        t(i)=d(i)
    end do
    if(err.lt.eps) exit
end do

do i=0,n
x=i*dx
    write(*,100) x, t(i), -1+sqrt(1+3.*x)
end do

! ファイルへの出力
open(1, file='Netsu.dat', status='replace')
do i=0,n
    x=i*dx
    write(1,100) x, t(i),-1+sqrt(1+3.*x)
end do
close(1)

100 format(2(f8.5,1x),f8.5)
end


!  TDMA
subroutine tdma(a,b,c,d,n)
real a(n), b(n), c(n),d(n)

! 前進消去
do i=1,n-1
    ratio=a(i+1)/b(i)
    b(i+1)=b(i+1)-ratio*c(i)
    d(i+1)=d(i+1)-ratio*d(i)
end do

! 後退代入
d(n)=d(n)/b(n)
do i=n-1,1,-1
    d(i)=(d(i)-c(i)*d(i+1))/b(i)
end do

end

 

なお、このプログラムでは、

  

と、点における熱伝導率の算術平均を用いて計算している。

熱伝導率のこの計算の是非については議論が必要だろうが、計算結果は極めて良好なのだから、「終わりよければ全てよし」ということで(^^)

 

本問は、いちおう、Tの非線形の微分方程式なので、計算においては、

  netsu-tei-0010.png

と、反復1回前の温度を用いて熱伝導率を計算し、線形化して解いている。

反復計算の収束条件は

  netsu-tei-0011.png

 


nice!(0)  コメント(0) 

問題(4月3日)の答らしきものとシュレーダー法 [数値解析]

問題(4月3日)の答らしきものとシュレーダー法

 

問題 次の方程式の解を2分法とニュートン法で求めたい。

  

さて、次の条件で計算したとき、2分法とニュートン法、どちらが速く収束解を得られるでしょうか。

 

条件

 2分法 計算開始の区間をI₀=[0,5]とする。

 ニュートン法 x₀=5を計算の出発点とする。

 収束判定 |x−3<10⁻⁶

 

【解答らしきもの】

  

したがって、

ニュートン法は

  

となる。

両辺から3を引くと

  

常用対数表.pngよって、収束条件を満たす自然数n

  

を満たさなければならない。

したがって、

  mon-0403-siki-001.png

ゆえに、n=36

 

一方、2分法の場合は、1回計算を進めるたびに方程式の解が存在する区間の幅が1/2になるので、

  

のとき、任意のに対して

  

となるので、、

  mon-0403-siki-002.png

最大、n=23回、計算すればよい。

 

故に、本問の場合、2分法の方がニュートン法よりも早く収束する。

(解答らしきもの終)

 

表計算ソフトで、ニュートン法を用いて解いた結果を以下に示す。

上の計算の通り、実際に、36回で収束条件を満たすことがわかるだろう。

 

次に、より一般の、x=αをm重根(m=2,3,・・・)とする

 

方程式の場合について考えることにする。

このとき、

  

になるので、ニュートン法は

  

となる。

この両辺をaで引くと、

  

したがって、m=2, 3,4,・・・と、mが大きくなればなるほど、ニュートン法は収束が遅くなることがわかる。

  

ここで、

  

を収束率と呼ぶことにする。

これは、1回計算すると、どれだけx=αという解に近づくいているかの比を表すもので、この値が小さければ小さいほど、速く収束することを表す。

(2)の場合、収束率は、厳密に

  

である。

(1)の場合、m=3だから、収束率は

  

となる。

 

さらに一般の、x=αm重根(m=2,3,・・・)とする次の方程式があるとする。

  

 

このとき、次のことが成り立つことが知られている。

 

定理

f(x)m重零点αを含む開区間Iにおいて級で、

  

とする。

にとり、

  

とすると、

  

である。

 

今回取り上げた問題は、その最も簡単なh(x)=1の場合。

 

そして、(4)式の方程式のm重根x=αを求めるとき、ニュートン法(5)を

  

と改良すると、2次の収束速度に早めることができることが知られている。

この方法をシュレーダー法という。

 

実際の方程式の近似解を求めるとき、シュレーダー法が役に立つかどうかというと、何とも微妙。何故ならば、求める解がm重根であるかどうか、我々は事前に知ることができないから(^^

m重根であるかどうかの判定をしながら計算する方法もあるらしいけれど、その判定、判断をするためにプログラムは複雑化し、その判定に要する時間も必要になる。結果、2分法や(5)式を用いたニュートン法よりも計算時間を要したのならば、いったい、何をしているのかがわからなくなってしまう!!

 

この他に、ハレー彗星で有名な、ハレー法(ハリー法)なども存在するようですが・・・。

 

 

定理の証明


nice!(0)  コメント(0) 

ニュートン法の収束速度 [数値解析]

ニュートン法の収束速度

 

§1 ニュートン法の計算法

 

x=αを解にもつ次の方程式があるとする。

  

f(x)が何回でも微分可能な関数であるとき、x=x₀f(α)をテーラー展開すると

  

となるので、2次の項を無視すると、f(x)を次のように近似できる。

  

x=αは(1)の解なのでf(α)=0であり、f(α)=0と(3)式より

  

と、(4)式を用いて方程式(1)の解を推測することができる。

(4)式の右辺で推測されたαの推測値をx₁とすると、

  

(5)式から得られた推測値x₁を用いて

  

と推測値x₂を求め、さらに、この推測値x₂をもとに・・・

  

と、逐次的に、方程式(1)の解の近似解を求める方法をニュートン法またはニュートン・ラプソン法という。

 

§2 ニュートン法の収束速度

 

f(α)でテーラー展開すると

  

ここで、ξαの間にある実数。

f(α)=0だから、のとき、

  

ここで、

   

とおくと、

  

となる。

ここで、適当なδ>0をとり、となるようにすると、仮定よりf'(x)f''(x)は、連続なので、Iで有界。

したがって、

  

となるA>0(仮定よりf'(x)>0)、B>0が存在する。

したがって、

  

となり、ニュートン法は2次の収束速度をもつことがわかる。

仮に、であるとすると、次のステップでは、その誤差の程度は0.01²=10⁻⁴になる。さらに、その次のステップでは、ニュートン法による誤差の程度は、(10⁻⁴)²=10⁻⁸と、急速に誤差が小さくなってゆく。

このため、ニュートン法は4、5回程度の反復計算で収束解を得ることができる。

 

 

f(x)=(x-3)^2(x+3).png上の表は、方程式

  

に対し、計算の初期値x₀=−5を与え、ニュートン法を用いて表計算ソフトで計算した結果である。5回の反復回数で、この方程式の解の1つであるx=−3が得られることがわかる。

また、反復計算の4、5回目の誤差の指数を見ると、−6と−12となっているとから、ニュートン法の収束速度が2次であることを確認できるだろう。

 

なのですが、

計算の初期値にx₀=5を与えると、事情は一変する。

 

表計算ソフトの誤差の列を見ると、反復計算の回数を1回増やすたびに誤差は約1/2にしかなっておらず、2次の収束速度を持っているはずのニュートン法が1次の収束速度に落ちてしまっていることがわかる。

このとき、ニュートン法は、2分法に限りなく近いものになってしまうのであった。

 

多重解x=αを持つ方程式

  

の解x=αをニュートン法で求めようとすると、一般に、非常に収束が遅い。

 

 

興味のあるヒトは、さらに、

  

の解x=3をニュートン法で求めてみるといいと思う。

 

 


nice!(0)  コメント(0) 
前の10件 | - 数値解析 ブログトップ