自分のブログログアウトねむねこ幻想郷ネムネコ・ファミリーのために、大昔に書いたファンタジー小説を公開。数学のサイトではありません。[ネコ騙し数学]ニュートン法による非線形連立方程式の近似解法2ニュートン法を用いると、非線形方程式の近似解(の一つ)をと求めることができる。修正するために必要なであるのだから、偏微分係数の値は近似値で十分なはずである。したがって、十分小さなhを選び、中心差分をもちいてと計算することにする。h=1/106程度にとれば、(6)式で計算した打切誤差は、せいぜい、1/1010程度で、その差はほとんど無視できる。丸め誤差の範囲の程度だにゃ。この方針に従ってあらたにプログラムを作ることにする。doend上のプログラムではとし、この連立方程式の解(の1つ)を求めている。計算開始の初期値(x0,y0,z0)=(1,1,1)として計算した結果は次の通り。[ネコ騙し数学]ニュートン法による非線形連立方程式の解法f(x,y)、g(x,y)を2回微分可能な関数とする。近似解法について考えることにする。(x,y)を連立方程式の解、(x0,y0)を(1)の解の推測値とする。f(x,y)、g(x,y)は仮定より微分可能なので、次のようにテーラー展開することが可能である。したがって、(x0,y0)が(x,y)の近傍の点とすると、と近似することができる。(x,y)は連立方程式(1)の解なので、f(x,y)=0、g(x,y)=0であるから、(4)、(5)式より、(1)式はと書き換えることが可能で、この連立方程式を解くことによって、(1)の解の推測値(x0,y0)をもとに解(x,y)の推測値を求めることができる。(6)式を行列を用いて書き換えると、ここで、とおくと、(8)の行列式|J|≠0のとき、(7)式よりと解くことができる。(9)式で得られた連立方程式の解(x,y)をあらたに(x0,y0)にし、(8)、(9)を再計算し、(x,y)を求める。この操作を繰り返すほど、(1)の近づくことが予想される。この方法は、そのまま、より一般のに拡張することが可能で、この場合は次のようになる。(12)から(非線形)連立方程式の近似解を求めることができるが、(12)を利用する場合、ヤコビ行列Jの逆行列J−1を求める必要がある。連立方程式(11)を解くよりヤコビ行列Jの逆行列を求めることの方が大変なので――連立方程式を直接解く場合と比較すると、Jの逆行列を求める計算は計算量が約n倍!!――なので、連立方程式をガウス消去法などで解き、と計算するとよい。このように(非線形)連立方程式の近似解を求める方法をニュートン(・ラプソン)法という。ニュートン法で解くプログラムを示す。u(x,t)=f(x−ct)であるとき、であることを示せ。ξ=x−ctとおき、合成関数の微分法を用いると、したがって、(解答終)したがって、fを微分可能な任意の関数とすると、u(x,t)=f(x−ct)はの(一般)解になる。f(x)=sinxとおくと、となり、これが(1)の(特殊)解になることから分かるように、(1)は波の方程式の一つである。(1)式の左辺第1項をと、時間微分に関しては前進差分を用いて近似し、(2)式の左辺の第2項をと、空間微分に関しては中心差分を用いて近似すると、(1)式はと近似することができ、これからという漸化式を得ることができる。(6)式の無次元数(物理の単位をもたない数のこと)をクーラン数という。いうことで、初期条件をさらに、c=0.2、Δx=1、Δt=1として、(1)を解いた結果は次のようになる。赤い曲線はx=1、黄色の曲線はx=2、緑はx=3、茶色はx=4におけるuの数値解の時間tにおける変化を表している。c=0.2だからx=1に高さ1の波が到達するのはt=5なのにも関わらず、t=Δt=1に既に波の一部(?)が到達しているだけでなく、t=5になっても波の高さは1にならない。波の高さは、約t=20以降、正弦曲線のような曲線を描き、振動する。x=2、x=3、x=4と下流にいけばいくほど、波の最大の高さは高くなり、より激しく振動するという、物理的にありえない、数値的に振動する解を(7)式は求めてしまう。Space)は、数値的に不安定な解法で、この解法を(1)式に用いると物理的にありえない振動解を求めてしまうということ。そこで、なぜ、FTCS法は不安定なのか、フォン・ノイマンの安定判定法を用いて調べることにする。し、(7)式に代入すると、となり、したがって、フォン・ノイマンの安定判定法によると、クーラン数に関わらず、FTCSは無条件不安定!!ということで、FTCSは微分方程式(1)の数値解を求めることに使えないことが明らかになった。そこで、(1)式の第2項の空間微分をと風上化(上流化)すると、c>0のとき、(1)はしたがって、となる。(8)式を(13)式に代入すると、よって、したがって、の時に安定になる。風上(上流)差分を用いて、c=0.2、Δx=1、Δt=1、すなわち、クーラン数C=0.2として計算した結果は、次の通り。高さ1の波の到達時刻は実際とは異なるけれど、今度は振動することなく計算できていることが分かる。下図を見ると、この波の解はでなければならないのに、風上差分を用いた(13)の数値解は、波の形が崩れていることが分かるだろう。参考までに、同一の条件でFTCSを用いて計算した結果を下図に示す。次の熱伝導方程式の右辺を差分法を用いて次のようと近似し、とおくと、ルンゲ・クッタ法を用いて、(1)式の近似解を求めることができる。例によって、(2)を用いて、次の問題を解くと、右のグラフのようになる。もとで、Δt=1、Δx=1として解け。Δt=1、Δx=1という粗い計算でも、良好な計算結果が得られていることが分かる。Δt=0.1、Δx=0.4とすると、さらに良好な結果が得られ、厳密解との差はほとんど認められなくなる。ルンゲ・クッタ法も陽解法の一種なので、のとき数値解が振動することがあるので、ルンゲ・クッタ法を用いて熱伝導方程式を解くとき、になるようにΔt、Δxを取ったほうがよい。このことを示すために、Δt=0.36、Δx=0.5のとき、すなわち、として計算すると、40ステップのt=14.4の前後で振動を始め、t=18ではさらに激しく振動する。計算に使用したプログラムは次の通り。[ネコ騙し数学][陽解法でも収束するかも]ddt^3です。面白いなぁ~と思って、差分法に関する安定条件をググってみると、フォン・ノイマンの安定条件ってのが出てきます(註1。陽解法が安定であるためには、任意の整数0≦mに対して、だそうです。ここでGは時間方向への誤差の拡大率,βmは、であり、Lは解析領域の長さです。この条件は解の大きさに誤差も比例するという前提なので、mは解を空間方向にフーリエ展開した場合の波数(モード次数)と考えられます。к=1/2,Δt=0.1,Δx=0.4,L=4として、m=0,1,2,・・・に対してG(m)をプロットしてみると、図-1になります。波数20以上は周期的な繰り返しです。解析結果のグラフをみてみると、どう考えても1次モードのみ励起されるはずです。Gは誤差の拡大率です。初期誤差をεとし、誤差は解の大きさに比例するという前提を信じれば、で、時間ステップnにおける誤差を評価できる事になります(たぶん最悪の予想)。T(x,t)は解,tnはnステップ目の時刻です。解析結果は放物線が一様に圧縮されてくような感じなので、T(x,tn)/T(x,0)をx=2の値の変化で代表させます。・・・表-1くらいかな?。これをグラフに起こして強引に指数関数近似したのが図-2です。相関係数がR2=1とウルトラ良いので(^^)、図-2の補間結果から、で近似する事にします(註2)。次に初期誤差εですが、T(x,0)=x(4-x)=-x2+4xなので、t=0でのx方向への2階差分の絶対値はa(0)=1のはず。時間差分は1階で、空間差分は2階だから、くらいかな?(^^;)。以上を使って、誤差の積算Ε(n)を計算してやると、図-3になります。陰解法は、長時間積分に対して陽解法より平均的に精度が良いだけの話で、局所的な精度は陽解法が上回る場合もあり得ると、自分は思っています。長くない時間積分であれば、無条件安定なシンプレクティック法よりルンゲクッタ法の方が、ずう~っとずう~っと精度が良い。シンプレクティック法は陽解法のくせに(まさに陰に)陰解法の特性を持つのに対して、ルンゲクッタ法による超長時間積分では、解は0に減衰します。予想1)陽解法でも発散しないケースだったから、陽解法の精度がバカみたいに良かった?。n=100で既に0.001~最大0.01程度の誤差は見込まれます。問題のスケール規模は100=1のオーダーで単精度の有効数字は6桁。倍精度でやりたいなぁ~(^^;)。予想2)倍精度計算したら、劇的に状況が変化するかも。以上、他人のプログラムなんて頼まれても読みたくないddt^3の、他人事のような発言でした(^^;)。(執筆:ddt3さん)(註1)フォン・ノイマン法による陽解法の安定性判定については、ねこ騙し数学でも過去に「陽解法による拡散方程式の安定性」という記事で取り上げている。http://nekodamashi-math.blog.so-net.ne.jp/2016-11-23(註2)微分方程式の厳密解は、そして、n=1のときとなるので、約5%くらいの誤差でと近似できる。まとめ熱伝導方程式を、陽解法を用いて離散化すると、(純)陰解法を用いると、となる。クランク・ニコルソン法の場合は、(2)、(3)式の辺々を足しあわせ2で割ったである。ここで、また、(2)、(3)、(4)式はと一つの式表すことができ、になる。クランク・ニコルソン法と解法の区別をすることなく、三重対角行列を係数に持つ連立方程式として(1)の近似解を求めることができる。この方針にしたがって作ったプログラムは以下の通り。クランク・ニコルソン法といった解法の区別をすることなく計算できる。αを一種の重みと考え、αを0、1/2、1だけではなく、0≦α≦1の任意のαを選び計算することもできる。とると、の値によって、伝播誤差のために、振動解が得られる。例によって、(1)のκ=1/2とし、条件とし、Δt=0.1、Δx=0.4に選んで計算した結果を以下に示す。クランク・ニコルソン法ともに、(1)の厳密解とよく一致することが分かる。Δt=0.2、Δx=0.4とし、陽解法で解いた結果は次の通り。このとき、となるので、この条件で陽解法を用いて解こうとすると、伝播誤差のために、計算を進めるたびに、数値解の振動が激しくなり、厳密解からの逸脱が大きくなることが分かる。微分方程式を、陽解法を用いて離散化すると、となる。陰解法を用いて離散化すると、(2)と(3)の平均を取ると、になる。ここで、とおくと、と3重対角行列を係数に持つ連立方程式が得られ、三重対角行列アルゴリズムTDMAを用いてこの連立方程式を解くことが可能になる。以外の項を落とすと、となり、r>0のとき、となり、クランク・ニコルソン法は無条件安定であることが予想できる。例によって、次の問題をクランク・ニコルソン法で解いてみる。もとで、Δt=1、Δx=1として解け。Δt=1、Δx=1という粗い計算でも、良好な計算結果が得られていることが分かる。Δt=0.5、Δx=0.4にすると、厳密解とクランク・ニコルソン法による近似解はほとんど等しくなる。この計算に使用したプログラムは次の通り。数値解法2前回に引き続き、差分法を用いた熱伝導方程式の数値解法について述べることにする。この方程式の時間微分を右辺の空間微分をと近似すると、(1)式はとなる。とおくと、という3重対角行列を係数を持つ連立方程式が得られる。(1)の解をT=T(x,t)とし、初期条件を境界条件をとする。すると、初期条件より境界条件よりとなる。したがって、という連立方程式が得られ、逐次的に、3重対角行列アルゴリズム(TDMA)を用いて解くことによって、微分方程式(1)の近似解を求めることができる。このような解法を陰解法という。陰解法は、前回紹介した陽解法とは異なり、無条件安定である。もとで、Δx=1、Δt=1として解け。解かせたものは次の通り。この問題の場合、厳密解と比較すると、陰解法による近似解は厳密解よりも大きいことがわかる。前回の陽解法による結果を加え比較したものが次のグラフである。込み入っていて、少しわかりづらいと思うのだけれど、厳密解は陰解法と陽解法の間に位置している。いうことで、陰解法と陽解法によって得られた数値解を足して2で割ると厳密解に近い結果が得られることが予想できる。この予想は正しく、陽解法と陰解法の中間的な解法、クランク・ニコルソン法という解法が存在する。Δt=1、Δx=1という非常に粗い計算でありながら、クランク・ニコルソン法は、厳密解とよく一致していることが分かる。この事実をもって、陽解法や陰解法は精度が悪いと思わないで欲しい。Δt=0.1程度にとれば、陽解法、陰解法ともに良好な結果を得ることができる。陽解法はともかく、陰解法は熱や流れなどの数値計算ではよく使われている、優れた解法である。問題の初期条件を境界条件をと変え、陰解法を用いてΔt=0.5、Δx=0.4として計算した結果を以下に示す。このようにあらい格子でも、厳密解をよく再現していることが分かる。t→∞の定常解に収束することが確認できる。誤差は最大で10%くらい。Δt=0.1、Δx=0.4にすると、誤差は最大で約3%。陰解法よりより高精度な計算ができる、クランク・ニコルソン法は、次回ということで。計算に使用したFortranプログラムは次の通り。熱伝導方程式を例に取り、差分法を用いた数値的な解法について考えることにする。右辺の2階の空間(偏)微分をと、差分法を用いて近似すると、ここで、とおくと、となる。計算することができる。r=1/2のとき、(5)式はになる。このような解法を陽解法という。陽解法は代数方程式を解くことなく簡単に計算できるが、のとき、(5)を用いて求めた近似解は安定ではなく、振動解が得られる。したがって、となるようにΔt、Δxを選ばないといけない。誤差がこれにしたがって伝播するとする、安定であるためには、でなければならい。したがって、安定であるためには・・・。すると、だから、少なくとも、安定であるためにはでなければならないに違いない。この条件を満たさないと、近似解は振動したり、発散するだろう。Neumannの判定法などを用いる必要があるが、それは厄介なので、簡易的にこの関係を求めてみたにゃ。もとで、Δx=1として解け。Δtをとなるように、Δt=1にとると、となる。以下、同様に計算すると、次の表が得られる。この結果をグラフにすると、次のようになる。定量的にはともかく、このような粗い計算であっても、指数関数的な現象を捉えており、定性的には正しい結果が得られていることがわかる。精度よく計算するために、Δx=0.1とすると、条件(7)よりよって、Δtを0.01以下にとる必要があり、計算量が大きく増えてしまう。したがって、実際は、熱伝導方程式を陽解法を用いて解くことはない。した制約のない陰解法を次回紹介することにする。doend出力部は自分で作るといいにゃ。そこまでは面倒を見きれないケロよ。mはX方向の分割数、nはY方向の分割数、n1=n+1だケロ。Xの最大値は10だから、m=100、n=10くらいが丁度いいんじゃないかな。mを1000、nを100にしても、m×nオーダーの計算法だから、エンターキーを押した瞬間に終わると思うけれど、U、V、T合せて30万のデータ数になるから、mとnをあまり大きくすのはやめたほうがいいと思うにゃ。今,確..コメント、ありがとうございます。..Cもお願いしますベートーヴェンについて。ベートーヴェン..ルートヴィヒ・ヴァン・ベートーヴェ..コメントありがとうございます。ソネットポイントこの広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。