[自動轉寄] Re: [程式] 有關mathematica求解的一個問題

作者: af2486 (我喜歡旅行)
標題: Re: [程式] 有關mathematica求解的一個問題
時間: Sat Jul 31 23:45:57 2010

您好~不好意思我想請問您的程式碼是要從哪邊開始貼呢?
或是說取代哪邊的程式碼...我取代"Do"這個指令後面的
要執行的時候好像有點小問題(出在函數的引數[]那邊,有個錯誤訊息)
不好意思才剛開始學連基本的debug都不太會...這小問題可能還要麻煩你了謝謝!

※ 引述《chungyuandye (養花種魚數月亮賞星星)》之銘言：
: ※ 引述《af2486 (我喜歡旅行)》之銘言：
: : ------------------------------------------------------------------------
: : [軟體程式類別]:Mathematica
: : 請填入軟體程式類別 例如SAS、SPSS、R、EVIEWS...等
: : [程式問題]:
: : 在用mathematica牛頓法求解時，會先給個固定值讓程式在那值附近求解(一連串的多組解)。
: : 不過後來發現往往後面求出的解會跑掉，因此現在要設定每次求解都以上一個算出來的值
: : 當作新的求解範圍。不過我不知道該用什麼指令好，有人能給個方向嗎?
: : xiLimit = 15.;
: : Hankel1[m_, x_] :=
: : If[Abs[x] < xiLimit, BesselJ[m, x] + I BesselY[m, x],
: : Sqrt[2/(Pi*x)] Exp[I*(x - m/2*Pi - Pi/4)]];
: : BesselJp[m_, x_] :=
: : If[Abs[x] < xiLimit,
: : D[BesselJ[m, x1],
: : x1] /. {x1 -> x}, -Sqrt[2/(Pi*x)] Sin[(x - m/2*Pi - Pi/4)]];
: : BesselYp[m_, x_] :=
: : If[Abs[x] < xiLimit, D[BesselY[m, x1], x1] /. {x1 -> x},
: : Sqrt[2/(Pi*x)] Cos[(x - m/2*Pi - Pi/4)]];
: : Hankel1p[m_, x_] :=
: : If[Abs[x] < xiLimit, D[Hankel1[m, x1], x1] /. {x1 -> x},
: : I*Sqrt[2/(Pi*x)] Exp[I*(x - m/2*Pi - Pi/4)]];
: : kr[kz_, k_] := Sqrt[k^2 - kz^2];
: : \[Alpha]r[kz_, k_] := Sqrt[kz^2 - k^2];
: : f[kz_, opt_] := Module[{sub},
: : k0 = \[Omega] Sqrt[\[Epsilon]0 \[Mu]0] /. opt;
: : k1 = \[Omega] Sqrt[\[Epsilon]1 \[Mu]1] /. opt;
: : k2 = \[Omega] Sqrt[\[Epsilon]2 \[Mu]2] /. opt;
: : r0 = r0 /. opt;
: : r1 = r1 /. opt;
: : m = m /. opt;
: : kr0 = kr[kz, k0];
: : kr1 = kr[kz, k1];
: : \[Alpha]r2 = \[Alpha]r[kz, k2];
: : BJ00 = BesselJ[m, kr[kz, k0] r0];
: : BJ10 = BesselJ[m, kr[kz, k1] r0];
: : BJ11 = BesselJ[m, kr[kz, k1] r1];
: : BY10 = BesselY[m, kr[kz, k1] r0];
: : BY11 = BesselY[m, kr[kz, k1] r1];
: : BJp00 = BesselJp[m, kr[kz, k0] r0];
: : BJp10 = BesselJp[m, kr[kz, k1] r0];
: : BJp11 = BesselJp[m, kr[kz, k1] r1];
: : BYp10 = BesselYp[m, kr[kz, k1] r0];
: : BYp11 = BesselYp[m, kr[kz, k1] r1];
: : HA21 = Hankel1[m, I \[Alpha]r[kz, k2] r1];
: : HAp21 = Hankel1p[m, I \[Alpha]r[kz, k2] r1];
: : M = {{BJ00, 0, -BJ10, -BY10 , 0, 0, 0, 0},
: : {0, BJ00, 0, 0, -BJ10, -BY10, 0, 0},
: : {(\[Omega] \[Epsilon]0)/(kr0 r0) BJp00, (m kz)/(kr0^2 r0^2)
: : BJ00, -((\[Omega] \[Epsilon]1)/(kr1 r0))
: : BJp10, -((\[Omega] \[Epsilon]1)/(kr1 r0)) BYp10, -((m kz)/(
: : kr1^2 r0^2)) BJ10, -((m kz)/(kr1^2 r0^2)) BY10, 0, 0},
: : {(m kz)/(kr0^2 r0^2) BJ00, (\[Omega] \[Mu]0)/(kr0 r0)
: : BJp00, -((m kz)/(kr1^2 r0^2)) BJ10, -((m kz)/(kr1^2 r0^2))
: : BY10, -((\[Omega] \[Mu]1)/(kr1 r0))
: : BJp10, -((\[Omega] \[Mu]1)/(kr1 r0)) BYp10, 0, 0 },
: : {0, 0, BJ11, BY11 , 0, 0, -HA21, 0},
: : {0, 0, 0, 0, BJ11, BY11, 0, -HA21},
: : {0, 0, (\[Omega] \[Epsilon]1)/(kr1 r1)
: : BJp11, (\[Omega] \[Epsilon]1)/(kr1 r1) BYp11, (m kz)/(
: : kr1^2 r1^2) BJ11, (m kz)/(kr1^2 r1^2) BY11, (
: : I \[Omega] \[Epsilon]2)/(\[Alpha]r2 r1) HAp21, (
: : m kz)/(\[Alpha]r2^2 r1^2) HA21},
: : {0, 0, (m kz)/(kr1^2 r1^2) BJ11, (m kz)/(kr1^2 r1^2)
: : BY11, (\[Omega] \[Mu]1)/(kr1 r1) BJp11, (\[Omega] \[Mu]1)/(
: : kr1 r1) BYp11, (m kz)/(\[Alpha]r2^2 r1^2) HA21, (
: : I \[Omega] \[Mu]2)/(\[Alpha]r2 r1) HAp21}
: : } /. opt;
: : equ = Det[M];
: : Return[equ];
: : ]
: : Do[
: : n0 = Sqrt[
: : 1 + (0.68671749* lambda^2)/(lambda^2 -
: : 0.072675189^2) + (0.43481505* lambda^2)/(lambda^2 -
: : 0.11514351^2) + (0.89656582* lambda^2)/(lambda^2 -
: : 10.002398^2)];
: : n1 = Sqrt[
: : 1 + (0.6961663* lambda^2)/(lambda^2 - 0.0684043^2) + (0.4079426*
: : lambda^2)/(lambda^2 - 0.1162414^2) + (0.8974794*
: : lambda^2)/(lambda^2 - 9.896161^2)];
: : n2 = 1.0;
: : opt = {m -> 1, \[Omega] -> 2 Pi/lambda, r0 -> 4.1,
: : r1 -> 62.5, \[Epsilon]0 -> n0^2, \[Mu]0 -> 1., \[Epsilon]1 ->
: : n1^2, \[Mu]1 -> 1., \[Epsilon]2 -> n2^2, \[Mu]2 -> 1., kz -> y};
: : sol2 = FindRoot[f[y*(2 Pi/lambda), opt] == 0, {y, 1.4435}];
: : Print[FullForm[Re[sol2[[1, 2]]]]],
: : {lambda, 1.54, 1.56, 0.001}]
: : 以上是部分程式碼，從1.54~1.56每次加0.001在1.4435附近求解，會求出21組解。現在想要
: : 每算完一組解就以這個解當新的值繼續往下算；例如說1.54算出1.4480196760169795，
: : 1.541就以這個值往下算，之後的以此類推。請問這樣程式碼要怎改比較好呢?謝謝!
: : [軟體熟悉度]:
: : 低(1~3個月)
: : -----------------------------------------------------------------------------
: n0[lambda_]:=Sqrt[1+(0.68671749*lambda^2)/(lambda^2-
: 0.072675189^2)+(0.43481505*lambda^2)/(lambda^2-
: 0.11514351^2)+(0.89656582*lambda^2)/(lambda^2-10.002398^2)];
: n1[lambda_]:=Sqrt[1+(0.6961663*lambda^2)/(lambda^2-0.0684043^2)+(0.4079426*
: lambda^2)/(lambda^2-0.1162414^2)+(0.8974794*
: lambda^2)/(lambda^2-9.896161^2)];
: n2[lambda_]:=1.0;
: opt[lambda_]:={m->1,\[Omega]->2Pi/lambda,r0->4.1,
: r1->62.5,\[Epsilon]0->n0[lambda]^2,\[Mu]0->
: 1.,\[Epsilon]1->n1[lambda]^2,\[Mu]1->1.,\[Epsilon]2->
: n2[lambda]^2,\[Mu]2->1.,kz->y};
: NestList[{#[[1]]+0.01,#[[3]],Re@y}/.
: FindRoot[f[y*(2Pi/(#[[1]]+0.01)),opt[#[[1]]+0.01]],{y,#[[3]]}]&,
: {1.53,1.4435,1.4435},20][[2;;-1]]//TableForm

--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 111.252.213.103