------------------------------------------------------------------------ [軟體程式類別]: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個月) ----------------------------------------------------------------------------- --

※ 發信站: 批踢踢實業坊(ptt.cc) ◆ From: 140.115.41.12 ※ 編輯: af2486 來自: 111.252.204.176 (07/30 23:00)