4変数の非線形方程式のときかた

質問

初めまして＾＾
現在大学の研究で素子の設計を行っているのですが、それに必要な計算式の解き方がわからずこまっています＞＜；
計算にはMatlabまたはScilabを用いてプログラムを書いています。
どなたか数値計算に強いかたがいればご教授願えないでしょうか？
よろしくお願いします。

数式とプログラムは長くなるため、省きます。

回答No.3

Kules · ベストアンサー

A No.1のKulesです。
>fsolveという関数があり、１変数の非線形方程式の解で、与えられた初期値に最も近い解を導出してくれます。
調べてみたらfsolveはoptimization toolboxを持ってないと使えないみたいですね。
私は持ってないので使えないです(笑)
１変数ならfzeroとかいう関数もあるみたいなんですけどね…

>ニュートン法は２次のテイラー展開を行って連立方程式を作り解く方法ですよね？
あれ…私の理解ではニュートン法は１次のテイラー展開
（って言うのかな？１次導関数までしか使わないやつです）
で片付けるやつだったと思うのですが…
↑あまりこの辺の難しい数学は知らないのです。
また、使うのも全微分ではなく偏微分なのでそんなに項が爆発的に増えることはないと
思うのですが。ああ、積の導関数だと積の数だけ増えますね
（f=suv(s,u,vはxの関数)となっていた場合、f'の項は３つになりますね）

ただ、ネットで検索していく中のどっかのページに微分出来ない関数
（微分がメンドクサイ関数でもいいと思います）を
df/dx=(f(x+h)-f(x))/h
で解くケースが見られましたので、偏微分に関しても先の例でいけば、
f1(a,b)をx1で微分したものにx1=a,x2=bを代入したものd/dx1(f1(a,b))を
d/dx1(f1(a,b))=(f1(a+h,b)-f1(a,b))/h
として解いても良さそうです。
ということは、f1,f2に当たる関数をMatlab内で作っておけば、ヤコビ行列
[df1/dx1,df1/dx2;df2/dx1,df2/dx2]
は
[(f1(a+h,b)-f1(a,b))/h,(f1(a,b+h)-f1(a,b))/h;(f2(a+h,b)-f2(a,b))/h,(f2(a,b+h)-f2(a,b))/h]
とできるので、そこまで複雑にはならなさそうですね。
hは1e-4でも何でも小さい数を適当に与えればよいと思います。

というかここまでニュートン法でつっぱらずに他の方法探れよって言われそうですが、
上にも書いたようにMatlabは触っているもののそこまで数学(特に大学以降の)には
詳しくないので…

参考になれば幸いです。

回答No.2

Tacosan

MATLAB は使ったことないので知りませんが, 検索したところによると
fsolve で多変数の方程式が解ける
ように読める.
万一 Optimization Toolbox がなくても, 検索すれば何とかなりそう.

参考URL：

http://www.mathworks.com/help/ja_JP/toolbox/optim/ug/fsolve.html

回答No.1

Kules

どんな数式なのかによってテクニックがあったりなかったりするかと思うんですが。
とりあえず泥臭い方法を書いてみます。
（Matlabにこの辺を扱える関数があったような気がしますが、私は使ったことないです）

ニュートン-ラフソン法を多変数に拡張してやれば解けるんじゃないかと思います。
１変数の時と変わるのは
・微分ではなく偏微分を使う
・割り算が行列の割り算（逆行列を掛ける）になる
ぐらいでしょうか。

４変数は書いてられないので、２変数で。
f1(x1,x2)=0、f2(x1,x2)=0を解くことを考えます。
適当な初期値a,bとして、それと解とのずれをΔx1、Δx2とします。
f1(a+Δx1,b+Δx2)=0、f2(a+Δx1,b+Δx2)=0
となればよいのですが、テイラー展開（だったかな？df1/dx1はf1をx1で微分することを表します）を
使って展開し、高次の微分の項を無視すれば
f1(a+Δx1,b+Δx2)=f1(a,b)+(df1/dx1)(Δx1)+(df1/dx2)(Δx2)
f2(a+Δx1,b+Δx2)=f2(a,b)+(df2/dx1)(Δx1)+(df2/dx2)(Δx2)
となります。
左辺=0として、右辺を行列の形で書くと上式は
[f1(a,b);f2(a,b)]+[df1/dx1,df1/dx2;df2/dx1,df2/dx2]*[Δx1;Δx2]=0
となるので、ここからΔ=[Δx1;Δx2]を求めることができます。
あとはa=a+Δx1、b=b+Δx2としてa,bの初期値を更新し、これを繰り返せば
初期値がいいところにあれば解は収束します。うまくいかなければ発散します。
変数が４つになってもやることは同じです。

Matlabでどう組むかは人によってそれぞれですね。ベタに書くならforループ１つで書けますが、
少しでも楽したいなら最低限f1,f2の関数は作った方がいいでしょう。

私なら関数f1(x1,x2),f2(x1,x2),df11(x1,x2),df12(x1,x2),df21(x1,x2),df22(x1,x2)
を作るかな。あるいはニュートン-ラフソン法以外の方法を探すかも知れません。
他になさそうならあきらめてこれでやるかと。

参考になれば幸いです。

お礼

丁寧な回答ありがとうございました。
ニュートン法自体はしっているのですが、4変数となると・・・「？？？」って感じです（汗）
ニュートン法は２次のテイラー展開を行って連立方程式を作り解く方法ですよね？
ニュートン法を4変数に拡張するためには、４変数のテイラー展開を求める必要があるかと思います。
頑張れば導出できるかもしれませんが（自信は全くないですが・・・；ω；）、複雑で長い計算式のためできれば微分をしない方法はないでしょうか？
この方法だと、100～200項くらいは出てきそうなので・・・・。

丁寧なご回答ありがとうございました。

補足

補足：
>-（Matlabにこの辺を扱える関数があったような気がしますが、私は使ったことないです）
fsolveという関数があり、１変数の非線形方程式の解で、与えられた初期値に最も近い解を導出してくれます。