• ベストアンサー

SOR法による差分方程式の解法

プログラミングのカテゴリーか少し悩んだのですが、 物理関係で差分方程式を解いている人が多いかと思い、 こちらで質問させて頂きます。 解きたいのは、2Dポアソン方程式です。(1Dでも構いません) 微分方程式を差分化してSOR法を使用します。 SOR法は若干ながらメッシュサイズ、加速係数に 解(ポテンシャル形状)が依存します。 しかし、できるだけこういった条件に依存せずに、 正確な解を得たいです。 より安定に正しい(メッシュサイズなどに依存しない)解が得られる 方法など、SOR法について何かご存知の方、何でも構いませんので、 ご教授のほど宜しくお願いします。 又、複雑な構造(ダブルHEMT構造、局所的に高ドープなど)でも、計算を 収束し易くする方法など、何かありましたらアドバイスをお願いします。

質問者が選んだベストアンサー

  • ベストアンサー
回答No.1

イオンビームの計算で2Dポアソン方程式をSOR法で解いたことがありますが、加速係数の決め方は試行錯誤でやってました。ポアソン方程式ですと右辺が問題依存なので教科書どおりにはいかず、加速係数を変えて反復回数が少なくなるようにするしかありません。面倒ですが、問題の規模が大きくなると反復回数が数千とか数万を超えるので加速係数を使うメリットはあります。 メッシュサイズについては、メッシュサイズを変えて解の収束性を調べ、自分がよしとする精度の解が得られるサイズを選ぶしかありません。どこまでをもってよしとするかはこれもまた問題依存ですが、計算機の性能の制約で決まる場合が多いです。 ところで、SOR法とメッシュサイズの問題は切り分けた方がいいでしょう。そもそも差分化によりポアソン方程式を有限な離散点で近似しているため、仮に差分式を何らかの方法で解いたとしてもそれがもとのポアソン方程式の正確な解であるかは別問題です。

gandhi-
質問者

お礼

貴重なご意見、有難うございます。大変参考になります。 >自分がよしとする精度の解が得られるサイズを選ぶしかありません。 そうですか、やはりバイアス、ドナー密度といった計算条件に合わせて、 適切に決めていくしかないようですね。 この設定で大丈夫、というメッシュの決め方があれば、相当プログラムが 使い易くなるので、もしあれば…と思ったんですが、そんな簡単にはいきませんか…。

全文を見る
すると、全ての回答が全文表示されます。

その他の回答 (1)

  • device
  • ベストアンサー率42% (3/7)
回答No.2

微分方程式を差分化して解いているということは、既にそこに近似的要因を含んでいることになります。よって、「正確な解」というものは得られません。「正確な解」を求める場合、解析的に解く必要があります。差分法で求める解は、「正確な解」と比較して、誤差の小さい解は求めることはできます。ご存知かとは思いますが、そのことには注意してください。 SOR法についてですが、加速係数は最適な値を求めることができます。このことに関しては、一般の文献にも書かれてあるのでそちらを参考にしてください。 メッシュサイズですが、場所によってメッシュサイズを変更するようなプログラミングをオススメします。変化の大きい部分に対しては、その変化に対して十分に小さいメッシュサイズを与えることによってうまく収束すると思います。この関数的なメッシュサイズの与え方によって、収束しなかったものも収束することがあります。 ポアソン方程式を解かれているようですが、その場合には、差分法によるプログラミングよりも、NEWTON法によるプログラミングを強くオススメします。計算時間が数万分の1になることも多々あります。

gandhi-
質問者

お礼

色々と教えてくださり、ありがとうございます。 >微分方程式を差分化して解いているということは、 >既にそこに近似的要因を含んでいることになります。 そうですね、その点は注意しておかないとダメですね。 >加速係数は最適な値を求めることができます。 図書館で何か本をあたってみようと思います。 >場所によってメッシュサイズを変更するような >プログラミングをオススメします。 はい、現在のプログラムは10レイヤー構造で、各レイヤーで メッシュサイズが設定できます。そうですね、ポテンシャル変化の 激しいレイヤーのメッシュを細かくすると収束し易いです。 回答頂きありがとうございましたm(__)m

全文を見る
すると、全ての回答が全文表示されます。

関連するQ&A

  • 差分法の精度(電位の計算)

    よろしくお願いします。 現在、差分法により電位を求めています。 ポアソン方程式を差分化してプログラムにより電位を求め、それと電位の式の計算結果をと比較しているのですが、差分法の計算結果があまり近い値となりません。 差分法の精度とはどのくらいなのでしょうか? 電荷がない場合では完璧に一致するのですが、電荷があると一致しません。ポアソン方程式の一般的な解はまだ導き出されておらず、電荷も常に一定ではないので十分な精度はない、みたいなことをちらっと聞いたこともあるので、もしやと思ったのですが、そもそも差分法の電位の計算とは一致しないものなのでしょうか? よろしくお願いします。

  • 連立方程式の解法

    連立方程式: Ax=B (A;係数行列、x;未知数、B;右辺行列) において detA = 0 であった場合、この解は一義的には定まらない という事なのですが、 このことはSOR法などの反復法も 使えないと言うこと言ってるのですか? detA = 0 の連立方程式はどうしても解けないのですか?

  • 偏微分方程式の数値解法

    偏微分方程式の込み入った質問です。 2次元(x,y)の空間で2つの関数f(x,y),g(x,y)を考えます。 そこで、それぞれにラプラス方程式を立てました。 fxx+fyy = 0  (1) gxx+gyy = 0 (2) です。これは境界値問題で、差分式からSOR法を使って収束計算によって数値解を求めることができます。f, gはそれぞれ独立という形にはなります。 そこにもう1つ式が出てきました。 fxfy + gxgy = 0 (3) というものです。f,gをx,yで1回微分してできる式です。 都合3つの式が出てきました。 この数値解を求めるにはどのような方法があるでしょうか。 数値解ですから近似解です。 3つ目の拘束条件の下でのラプラス方程式とみると、ペナルティ関数とかラグランジュの未定係数法とかいろいろあるかもなと思いますが。 3つ目の式は完全に満たすというより、できるだけ満足するようにしたいというものです。 よろしくお願いします。

  • ラプラス方程式の解法(情報処理として)と共役勾配法

    情報処理の数学はカテゴリ違いかもしれませんが、適当なところを見出せませんでしたので。 (以下、長文で申し訳ありません。) ラプラス方程式やポアソン方程式の境界値問題を解く方法にSOR法が あります。これは連立一次方程式の反復解法の1つだとも言えると思い ます。連立一次方程式の解法は通常AX=BとなってマトリックスAと既知 ベクトルBに対して未知ベクトルXを解くわけですが、ラプラス方程式 をSOR法で解く場合、AX=Bという形にはしません(そういう風に表現 はできますが)。例えば2次元を考えると100×100の2次元の領域を解く 場合、未知数は10000個です。これをAX=Bという風に考えたら、10000元 連立一次方程式を解くことになり、A(10000,10000)のようにメモリ確保 も大変ですが、実際にはそのように解きませんね。マトリックスAとい うメモリの取り方をしないのです。 100×100の格子領域でSOR法で解くとプログラムステップとしては数十 行ぐらいとなり、大したことはありませんね。 しかし、情報処理の連立一次方程式の解法の本を見ると、それはあくま でも連立1次方程式があった場合ということを前提としているのでその ようなメモリの確保から解説がスタートするようです。 そこで質問ですが、ラプラスやポアソン方程式を(前処理付)共役勾配 法などで解く場合、メモリ確保の前提条件としては連立一次方程式の形 ではないはずです。本を見るとマトリックス解法としての共役勾配法が 出ていますが、ラプラス方程式、ポアソン方程式の特化して解説してい る本などないでしょうか。また、簡単に説明できるのであれば教えて頂 きたいのですが。

  • 数値解析の手法(差分法)について

    現在、とある2元の1階偏微分方程式(解はu,vでそれぞれ右と左に進行する波)を数値解析によって解こうと考えています。 数値解析の手段として、差分法がよく用いられると思いますが、 現在、私は、場所に関してはuは後退差分、vは前進差分を使い、 時間に関しては前進差分を使って解いています。 ネットでは場所に関しては中心差分、時間に関してはルンゲ=クッタやリープフロッグなどが 使われていることが多く、私もこの2つを用いて解いてみました。 偏微分方程式には線形項が含まれていたため、 線形問題に対する制約であるΔt/Δx << 1は最低満たすように刻み幅をいろいろ取り、 計算時間も辞さず計算機を動かしてみましたが、 ノイズが消えず、解析解に限りなく近づくには至りませんでした。。 Δt/Δx=0.0001なども試したのですが・・・ そこで、諦めて違う差分法を試し、 場所に関して、uは後退差分、vは前進差分を使い、 時間に関しては前進差分を使って見たところ、 Δt/Δx=0.01程度で解析解に近い、なかなか精度の良い数値解を得ることが出来ました。 2次の差分では上手くいかず、1次の差分だとわりかし上手くいく・・・ 精度的には中心差分やルンゲ=クッタなどの方がいいと思うのですが・・ 正直不思議でなりませんでした。。。 最初に試した差分法のコードミスかと思い、何回もコードを確認し直しましたが、 やはり解析解に近づくには至りませんでした。 こんなことってあるのでしょうか?? 差分法でも場合によって使い分ける必要があるということでしょうか・・? その場合分けするときの指標など、知っておられる方、教えて頂けると助かります。 問題によって

  • 差分法でのメッシュ分割

    差分法を使用するときはまず領域(2次元)をメッシュ分割します。 C言語でプログラムを行う際に、メッシュ分割はどうプログラムしたらいいのでしょうか?今は二次元の配列を使用しています。 下図:二次元の配列a[3][4]の場合   [0][1][2][3] [0] □ □ □ □ [1] □ □ □ □ [2] □ □ □ □ これは、メッシュ分割したという風にいえるのでしょうか? 差分法でラプラス方程式を解いています。また、等間隔の格子間隔hを使用しています。 v(x+h,y)+v(x-h,y)+v(x,y+h)+v(x,y-h)-4v(x,y)=0 の連立方程式を解いています。

  • ナビィヴェ・ストークスの方程式の解法

    ナビエ・ストークスの方程式って解けるものなんですか? 差分法とか使わずに、単振動の線形微分方程式を解くかのごとく、手計算で厳密解を出したいのですが・・・(もちろん非線形微分方程式ですから線形と対比する事自体まちがっているかもしれませんが) 既存の研究ではレイノルズ数の小さいときには、定常解、非定常解ともに解の存在と一意性が証明されているらしいですが、その解を数式で記述することは可能なのでしょうか? 今は解けなくても、数十年後(未来)には解かれているのでしょうか? もしくは既存の関数だけでの記述は不可能で、何か今までに無い新しい関数を作らないと記述できないのでしょうか?(よく知らないけど、たしかディラックのデルタ関数なんかは出た当時は革命的なものだったでしょう) また、これらに関する最近の研究結果、論文等を知っている方、教えてください。 WEBのページ(日本語、英文以外はちょっと・・・)の紹介でもいいのですが、そのページに何が書かれていて、何が参考になるかなど書いてもらえると助かります。漠然とページを紹介されても、私は無知なんでそのページが何なのか理解できない事があります。  あと、質問内容自体、専門の方から見れば用語の取り違いとかあるかもしれませんが、質問の意味は伝わってますよね。私は一般人なのでそのへんは勘弁ください。

  • ラプラス方程式の解析解

    電磁気学を勉強しているのですが,分からないことがあるので質問させてもらいます. 静電場内にある電荷が作る電位分布を示す方程式としてラプラス方程式(∇^2*V=0)があると思います. ラプラス方程式とポアソン方程式の違いまでは理解できていると思います. 2次元のラプラス方程式は以下の式を変数分離法を用いて解くことで,直交座標系や球面座標系として考えることで,解析解が得られると理解しています. (ここまではたどり着くことが出来ました) (∂^2/∂^2x)V(x,y)+(∂^2/∂^2y)V(x,y)=0 分からないのは,ここから実際の電位分布を求める方法です. 具体的には,xy平面上の原点にポテンシャルV0がある場合,このV0による電位分布を求めることが出来ません. 直交座標系で考えると一般解は,A,B,C,D,kを定数として,次のようになると思います. V(x,y)=(A*exp(kx)+B*exp(-kx))*(Csinky+Dcosky) 境界条件から未知定数を求めたいのですが,うまくいきません・・・. 原点にポテンシャルがあるので,x→∞でV→0,y→∞でV→0,x=0,y=0でV=V0が境界条件になると思ったのですが,y→∞で(Csinky+Dcosky)は0に収束しません. 境界条件の設定が間違っているのでしょうか? 数値解では原点にポテンシャルを設定している解説は見つけられたのですが,解析解では資料がなく,どうすればいいか困っています. すみませんが,教えてください.

  • 拡散方程式の数値計算(有限差分法)による解き方

    下記の拡散方程式を数値計算(有限差分法)で解き方について教えてください。 拡散方程式 D∂^2C(x,y,t)/∂x^2+u∂C(x,y,t)/∂x=∂C(x,y,t)/∂t 境界条件 D∂C(x,y,t)/∂x=(K/β-u)×C(x,y,t)-KC'(y,t) at x=0 初期条件 C(x,y,0)=Co at t=0 質量保存則 ∂C'(y,t)/∂y=4/vd(K/β×C(0,y,t)-KC'(y,t) ) ---------------------------------------------------------- また、可能であれば、有限差分法以外にも数値解法(有限要素法、有限体積法など)があると思いますが、拡散方程式を解く場合、どの解法が一般的に適しているのか教えてください! よろしくお願い申し上げます。

  • 波動方程式の差分法による境界条件

    波動方程式(ρ*∂^2u/∂t^2=T*∂u^2/∂x^2)を差分法で、境界条件をx=0とx=Lで自由条件(T*∂u/∂x=0)とした場合を考えています。 上の波動方程式を差分化すると、 (u[n+1][j]-2*u[n][j]+u[n-1][j])/(dt^2)=(u[n][j+1]-2*u[n][j]+u[n][j-1])/(dx^2) の形になると思います。(T/ρ=1.0とし、nは時間、jは距離の格子点として考えています) 初期条件は適当な形状を与えます。 両端を固定条件(u[n][0]=0,u[n][xp]=0,)とした場合はうまく解を得ることが出来ました。(xpはx=Lでの格子点) 問題は自由条件(T*∂u/∂x=0)すなわち、 T*(u[n][1]-u[n][0])/dx=0、T*(u[n][xp+1]-u[n][xp])/dx=0 となる場合、これをどのように使用したらよいのでしょうか? または根本的に考え方が間違っているのでしょうか? 本当に困ってます。よろしくお願いいたします。 内容が不十分の場合は補足要求お願いします。