2010年7月31日 星期六

[自動轉寄] 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

2010年7月30日 星期五

Fwd: [程式] 有關mathematica求解的一個問題

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]]==
0,{y,#[[3]]}]&,{1.53,1.4435,1.4435},
20][[2;;-1]]//TableForm

---------- 轉寄的郵件 ----------
寄件者: chungyuandye.bbs@ptt.cc <chungyuandye.bbs@ptt.cc>
日期: 2010年7月31日上午2:55
主旨: [程式] 有關mathematica求解的一個問題
收件者: chungyuandye@gmail.com


作者: af2486 (我喜歡旅行) 看板: Statistics
標題: [程式] 有關mathematica求解的一個問題
時間: Thu Jul 29 17:26:11 2010


------------------------------------------------------------------------

[軟體程式類別]:Mathematica
 [1;30m請填入軟體程式類別  例如SAS、SPSS、R、EVIEWS...等 [m

[程式問題]:
在用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
 [1;37m推  [33mchungyuandye [m [33m:沒有f的定義                                        [m 07/30 10:54
※ 編輯: af2486          來自: 111.252.204.176      (07/30 23:00)
 [1;31m→  [33maf2486 [m [33m:已附上完整程式碼                                         [m 07/30 23:02




--
Chung-Yuan Dye

養花種魚數月亮賞星星
http://cydye1069.blogspot.com

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

作者: af2486 (我喜歡旅行) 看板: Statistics
標題: [程式] 有關mathematica求解的一個問題
時間: Thu Jul 29 17:26:11 2010


------------------------------------------------------------------------

[軟體程式類別]: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
推 chungyuandye:沒有f的定義  07/30 10:54
※ 編輯: af2486 來自: 111.252.204.176 (07/30 23:00)
→ af2486:已附上完整程式碼  07/30 23:02

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

作者: af2486 (我喜歡旅行) 看板: Statistics
標題: [程式] 有關mathematica求解的一個問題
時間: Thu Jul 29 17:26:11 2010


------------------------------------------------------------------------

[軟體程式類別]: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
推 chungyuandye:沒有f的定義  07/30 10:54
※ 編輯: af2486 來自: 111.252.204.176 (07/30 23:00)
→ af2486:已附上完整程式碼  07/30 23:02

2010年7月29日 星期四

Re: [問題] mathematica 的block用法和三個問題

作者: Passions (passion)
標題: Re: [問題] mathematica 的block用法和三個問題
時間: Sun Jul 25 16:39:26 2010

※ 引述《chungyuandye (養花種魚數月亮賞星星)》之銘言:
: ※ 引述《Passions (passion)》之銘言:
: : 謝謝你的回信!
: : 不過仍然有幾個問題:
: : 一、關於 Block 的用法,仍有一點不清楚。
: : 原式中的 Block[{A},B ; C ; D ]
: : 分號; 在 mathematica 中,應該是代表指令結束,
: : 但在這邊還是被包在 Block 裡面,這真的是非常奇怪,不知道語法是什麼。
: 忘了說,;在mathematica表示指令結束但不輸出,
: 不加分號就是指令結束而且輸出!
: qq=(a+b)^2;ss=(a+b)^2
: qq
: ss
: 但這個在Block、Modual、With裡面就不行
: 指令動作之前一定要加分號
: 有問題在到我blog留言討論吧∼

真的是太謝謝你了! 我好像看懂了!

若程式碼為: f[x_]:=Block[{A},B ; C ; D ]

則 Block 的 "body" 是 B;C;D (還是只有 B?不過感覺是BCD都屬於Block的body)

而 {A} 是用來形容哪些變數是不受外界影響的。


然後,因為只輸出D的值,所以也就等同於 f[x_]=D ,而 B、C 只是用來幫助運算

不曉得這樣子理解有沒有錯…


--
我有到你的部落格留言喔~ :)

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

(Passions) Re: [問題] mathematica 的block用法和三個問題

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [問題] mathematica 的block用法和三個問題
時間: Sat Jul 24 19:32:20 2010


※ 引述《Passions (passion)》之銘言:
: ※ 引述《chungyuandye (養花種魚數月亮賞星星)》之銘言:
: : {r,n,t,y}是區域變數,只有在f這個函數才有作用,一旦跳出f它什麼都不是
: : 舉個例子
: : g[x_]:=Block[{w=10,xx=x},2xx+w]
: : 執行一下
: : g[10]
: : {ww,xx}
: : g[10]=30
: : {ww,xx}={ww,xx}
: : 不過為了怕之前L,G,K這些變數已經有定義過,建議改成這樣
: : f[L_?NumberQ,G_?NumberQ,K_NumberQ]:=Block[{LL=L,GG=G,KK=K,r,n,t,y},
: : {r}=NDSolve[{包含n,t,LL,GG的微分方程式,n[0]==K},n,{t,0,10}];
: : y=n[a] /. r;Plus @@ ((b-y)^2)]
: : @@=Apply,Plus @@ ((b-y)^2)] =>把加法套用在((b-y)^2)上
: 謝謝你的回信!
: 不過仍然有幾個問題:
: 一、關於 Block 的用法,仍有一點不清楚。
: 原式中的 Block[{A},B ; C ; D ]
: 分號; 在 mathematica 中,應該是代表指令結束,
: 但在這邊還是被包在 Block 裡面,這真的是非常奇怪,不知道語法是什麼。

忘了說,;在mathematica表示指令結束但不輸出,
不加分號就是指令結束而且輸出!
qq=(a+b)^2;ss=(a+b)^2
qq
ss

但這個在Block、Modual、With裡面就不行
指令動作之前一定要加分號

有問題在到我blog留言討論吧∼


--
養花種魚數月亮賞星星

http://cydye1069.blogspot.com

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

(Passions) Re: [問題] mathematica 的block用法和三個問題

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [問題] mathematica 的block用法和三個問題
時間: Sat Jul 24 19:17:36 2010

※ 引述《Passions (passion)》之銘言:
: ※ 引述《chungyuandye (養花種魚數月亮賞星星)》之銘言:
: : {r,n,t,y}是區域變數,只有在f這個函數才有作用,一旦跳出f它什麼都不是
: : 舉個例子
: : g[x_]:=Block[{w=10,xx=x},2xx+w]
: : 執行一下
: : g[10]
: : {ww,xx}
: : g[10]=30
: : {ww,xx}={ww,xx}
: : 不過為了怕之前L,G,K這些變數已經有定義過,建議改成這樣
: : f[L_?NumberQ,G_?NumberQ,K_NumberQ]:=Block[{LL=L,GG=G,KK=K,r,n,t,y},
: : {r}=NDSolve[{包含n,t,LL,GG的微分方程式,n[0]==K},n,{t,0,10}];
: : y=n[a] /. r;Plus @@ ((b-y)^2)]
: : @@=Apply,Plus @@ ((b-y)^2)] =>把加法套用在((b-y)^2)上
: 謝謝你的回信!
: 不過仍然有幾個問題:
: 一、關於 Block 的用法,仍有一點不清楚。
: 原式中的 Block[{A},B ; C ; D ]
: 分號; 在 mathematica 中,應該是代表指令結束,

執行B->執行C->最後輸出D
你就當成你在寫程式就好了∼∼

: 但在這邊還是被包在 Block 裡面,這真的是非常奇怪,不知道語法是什麼。
: 二、Plus @@ ((b-y)^2)]
: 這個跑出來的結果會是 2 + b - y
: 這真的是非常奇怪,直接把它當2看,這樣不就把平方的意義給弄掉了嗎

這個一點都不奇怪,因為在Mathematica裡面
如果你把((b-y)^2)當成一個List,那第一個元素是b,第二個元素是-y,
第三個元素2,你把加法套到這三個元素當然是2+b-y

: 不過又回到第一個問題來,還是他在Block 裡能表現出特別的意義來?
: 三、舉的這個例子我能看懂
: 跑 g[x_]:=Block[{w=10,xx=x},2xx+w]
: g[10] 會得到30。
: 但是例子中的 {ww,xx}和 {ww,xx}={ww,xx}
: : g[10]
: : {ww,xx}
: : g[10]=30
: : {ww,xx}={ww,xx}
: 是什麼意思呢?
ww,xx都是區域變數,所以跳出g以外那根本沒有被定義
所以還是只傳回ww,xx

: 不好意思問題蠻多的 ><
: 如果能抽空幫忙回答的話,我會非常感激的
: 感謝!


--
養花種魚數月亮賞星星

http://cydye1069.blogspot.com

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

(Passions) Re: [問題] mathematica 的block用法和三個問題

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [問題] mathematica 的block用法和三個問題
時間: Sat Jul 24 07:15:47 2010

※ 引述《Passions (passion)》之銘言:
: 想請問 Block 的用法是什麼?
: 這一段程式碼我看了好久,也有找 help 來看,但一直都看不懂
: 可以請高手幫我解答一下嗎:
: f[L_?NumberQ,G_?NumberQ,K_NumberQ]:=Block[{r,n,t,y},
: {r}=NDSolve[{包含n,t,L,G的微分方程式,n[0]==K},n,{t,0,10}];
: y=n[a] /. r;Plus @@ ((b-y)^2)]
: ...
: 最後目的是對 f 這個自訂函數找 FindMinimum.
: 我的問題是:
: 一、Block 是不是用來創造出一個環境是使得裡面的數值不受外界影響,
: 然後指定不受影響的數為 r,n,t,y。不太能了解用在這邊的意義為何?
: 二、思考流程是什麼?是先從 NDSolve解出數值,再丟給 Block
: 做不知道什麼的用途,之後再傳給 f 函數?@@
: 三、Plus @@ 這個用法是怎麼去做計算的呢,我開新頁去測試,仍然理不出頭緒
: 請問有沒有大大能夠幫忙回答的呢?因為還蠻急的,找了不少人問都沒有結果。
: 感謝!

{r,n,t,y}是區域變數,只有在f這個函數才有作用,一旦跳出f它什麼都不是
舉個例子
g[x_]:=Block[{w=10,xx=x},2xx+w]
執行一下

g[10]
{ww,xx}

g[10]=30
{ww,xx}={ww,xx}

不過為了怕之前L,G,K這些變數已經有定義過,建議改成這樣
f[L_?NumberQ,G_?NumberQ,K_NumberQ]:=Block[{LL=L,GG=G,KK=K,r,n,t,y},
{r}=NDSolve[{包含n,t,LL,GG的微分方程式,n[0]==K},n,{t,0,10}];
y=n[a] /. r;Plus @@ ((b-y)^2)]

@@=Apply,Plus @@ ((b-y)^2)] =>把加法套用在((b-y)^2)上


--
養花種魚數月亮賞星星

http://cydye1069.blogspot.com

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

(af2486) Re: [程式] 有關mathematica求解的一個問題

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [程式] 有關mathematica求解的一個問題
時間: Fri Jul 30 08:56:42 2010

※ 引述《af2486 (我喜歡旅行)》之銘言:
: ------------------------------------------------------------------------
: [軟體程式類別]:Mathematica
: 請填入軟體程式類別 例如SAS、SPSS、R、EVIEWS...等
: [程式問題]:
: 在用mathematica牛頓法求解時,會先給個固定值讓程式在那值附近求解(一連串的多組解)。
: 不過後來發現往往後面求出的解會跑掉,因此現在要設定每次求解都以上一個算出來的值
: 當作新的求解範圍。不過我不知道該用什麼指令好,有人能給個方向嗎?
: 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個月)
: -----------------------------------------------------------------------------


沒有你的f,隨便假設一下

{#, y} /. FindRoot[n2 y^2 + 4 n1 y + n0 /. lambda -> #, {y, 1}] & /@
Range[1.54, 1.56, 0.0001]

--
養花種魚數月亮賞星星

http://cydye1069.blogspot.com

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

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

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [程式] 有關mathematica求解的一個問題
時間: Fri Jul 30 08:56:42 2010

※ 引述《af2486 (我喜歡旅行)》之銘言:
: ------------------------------------------------------------------------
: [軟體程式類別]:Mathematica
: 請填入軟體程式類別 例如SAS、SPSS、R、EVIEWS...等
: [程式問題]:
: 在用mathematica牛頓法求解時,會先給個固定值讓程式在那值附近求解(一連串的多組解)。
: 不過後來發現往往後面求出的解會跑掉,因此現在要設定每次求解都以上一個算出來的值
: 當作新的求解範圍。不過我不知道該用什麼指令好,有人能給個方向嗎?
: 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個月)
: -----------------------------------------------------------------------------


沒有你的f,隨便假設一下

{#, y} /. FindRoot[n2 y^2 + 4 n1 y + n0 /. lambda -> #, {y, 1}] & /@
Range[1.54, 1.56, 0.0001]

--
養花種魚數月亮賞星星

http://cydye1069.blogspot.com

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

2010年7月27日 星期二

Re: [問題] 請問插入eps但漏字的情形??

作者: chungyuandye (養花種魚數月亮賞星星) 看板: LaTeX
標題: Re: [問題] 請問插入eps但漏字的情形??
時間: Tue Jul 27 14:47:58 2010

※ 引述《roloc (回到原點)》之銘言:
: ※ 引述《marblepillar (櫻吹雪遍地殘紅)》之銘言:
: : 首先我用Mathematica畫圖存成eps檔,
: : 然後在Latex插入時發現有漏字的狀況@@
: : 舉例來說:
: : 圖中有五條線是不同參數q的結果,
: : 所以我圖中有標示q=1,q=2,q=3,q=4,q=5 對應到不同顏色的線,
: : 但是最後pdf只顯示q= ,q=2,q=3,q=4,q=5
: : 原本是1的地方變成空白了!
: 基本上建議使用 psfrag 這個 package,
: 也就是說,你在 mathematica 上註明 A B C D E,
: 然後在 tex 上下指令 \psfrag{A}{$q=1$} 等等。
: 如果還不太懂的話我再寫清楚一點。
: : 我試過把Mathematica畫的圖直接存成jpg檔就無此問題,
: : 但是jpg的解析度不夠好,所以還是想用eps。
: : 請問有辦法解決嗎?
: : 非常感謝囉~<(_ _)>

http://wwwth.mppmu.mpg.de/members/jgrosse/mathpsfrag/

有人在mathematica寫了一個的配合tex上psfrag的package
這個package會自動幫你把圖檔上文字產生psfrag的格式
在tex上直接include進去即可!

如果想要簡單一點的方法,在Mathematica中程式就先不要標示文字
程式跑完後再用drawing tools自己畫,也可以避免文字被檔到!

--
我打研究室走過 那獨坐電腦前的容顏如苦瓜的糾結
靈感不來 長壽的煙霧不散
研究室如小小的寂寞的城 恰如商管的電梯向晚

http://chungyuandye.twbbs.org

--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 218.173.131.13
推 marblepillar:感謝分享兩個方法:)  07/28 03:15

2010年7月25日 星期日

[自動轉寄] Re: [問題] mathematica 的block用法和三個問題

作者: Passions (passion)
標題: Re: [問題] mathematica 的block用法和三個問題
時間: Sun Jul 25 16:39:26 2010

※ 引述《chungyuandye (養花種魚數月亮賞星星)》之銘言:
: ※ 引述《Passions (passion)》之銘言:
: : 謝謝你的回信!
: : 不過仍然有幾個問題:
: : 一、關於 Block 的用法,仍有一點不清楚。
: : 原式中的 Block[{A},B ; C ; D ]
: : 分號; 在 mathematica 中,應該是代表指令結束,
: : 但在這邊還是被包在 Block 裡面,這真的是非常奇怪,不知道語法是什麼。
: 忘了說,;在mathematica表示指令結束但不輸出,
: 不加分號就是指令結束而且輸出!
: qq=(a+b)^2;ss=(a+b)^2
: qq
: ss
: 但這個在Block、Modual、With裡面就不行
: 指令動作之前一定要加分號
: 有問題在到我blog留言討論吧∼

真的是太謝謝你了! 我好像看懂了!

若程式碼為: f[x_]:=Block[{A},B ; C ; D ]

則 Block 的 "body" 是 B;C;D (還是只有 B?不過感覺是BCD都屬於Block的body)

而 {A} 是用來形容哪些變數是不受外界影響的。


然後,因為只輸出D的值,所以也就等同於 f[x_]=D ,而 B、C 只是用來幫助運算

不曉得這樣子理解有沒有錯…


--
我有到你的部落格留言喔~ :)

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

2010年7月24日 星期六

[自動轉寄] (Passions) Re: [問題] mathematica 的block用法和三個問題

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [問題] mathematica 的block用法和三個問題
時間: Sat Jul 24 19:32:20 2010


※ 引述《Passions (passion)》之銘言:
: ※ 引述《chungyuandye (養花種魚數月亮賞星星)》之銘言:
: : {r,n,t,y}是區域變數,只有在f這個函數才有作用,一旦跳出f它什麼都不是
: : 舉個例子
: : g[x_]:=Block[{w=10,xx=x},2xx+w]
: : 執行一下
: : g[10]
: : {ww,xx}
: : g[10]=30
: : {ww,xx}={ww,xx}
: : 不過為了怕之前L,G,K這些變數已經有定義過,建議改成這樣
: : f[L_?NumberQ,G_?NumberQ,K_NumberQ]:=Block[{LL=L,GG=G,KK=K,r,n,t,y},
: : {r}=NDSolve[{包含n,t,LL,GG的微分方程式,n[0]==K},n,{t,0,10}];
: : y=n[a] /. r;Plus @@ ((b-y)^2)]
: : @@=Apply,Plus @@ ((b-y)^2)] =>把加法套用在((b-y)^2)上
: 謝謝你的回信!
: 不過仍然有幾個問題:
: 一、關於 Block 的用法,仍有一點不清楚。
: 原式中的 Block[{A},B ; C ; D ]
: 分號; 在 mathematica 中,應該是代表指令結束,
: 但在這邊還是被包在 Block 裡面,這真的是非常奇怪,不知道語法是什麼。

忘了說,;在mathematica表示指令結束但不輸出,
不加分號就是指令結束而且輸出!
qq=(a+b)^2;ss=(a+b)^2
qq
ss

但這個在Block、Modual、With裡面就不行
指令動作之前一定要加分號

有問題在到我blog留言討論吧∼


--
養花種魚數月亮賞星星

http://cydye1069.blogspot.com

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

[自動轉寄] (Passions) Re: [問題] mathematica 的block用法和三個問題

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [問題] mathematica 的block用法和三個問題
時間: Sat Jul 24 19:17:36 2010

※ 引述《Passions (passion)》之銘言:
: ※ 引述《chungyuandye (養花種魚數月亮賞星星)》之銘言:
: : {r,n,t,y}是區域變數,只有在f這個函數才有作用,一旦跳出f它什麼都不是
: : 舉個例子
: : g[x_]:=Block[{w=10,xx=x},2xx+w]
: : 執行一下
: : g[10]
: : {ww,xx}
: : g[10]=30
: : {ww,xx}={ww,xx}
: : 不過為了怕之前L,G,K這些變數已經有定義過,建議改成這樣
: : f[L_?NumberQ,G_?NumberQ,K_NumberQ]:=Block[{LL=L,GG=G,KK=K,r,n,t,y},
: : {r}=NDSolve[{包含n,t,LL,GG的微分方程式,n[0]==K},n,{t,0,10}];
: : y=n[a] /. r;Plus @@ ((b-y)^2)]
: : @@=Apply,Plus @@ ((b-y)^2)] =>把加法套用在((b-y)^2)上
: 謝謝你的回信!
: 不過仍然有幾個問題:
: 一、關於 Block 的用法,仍有一點不清楚。
: 原式中的 Block[{A},B ; C ; D ]
: 分號; 在 mathematica 中,應該是代表指令結束,

執行B->執行C->最後輸出D
你就當成你在寫程式就好了∼∼

: 但在這邊還是被包在 Block 裡面,這真的是非常奇怪,不知道語法是什麼。
: 二、Plus @@ ((b-y)^2)]
: 這個跑出來的結果會是 2 + b - y
: 這真的是非常奇怪,直接把它當2看,這樣不就把平方的意義給弄掉了嗎

這個一點都不奇怪,因為在Mathematica裡面
如果你把((b-y)^2)當成一個List,那第一個元素是b,第二個元素是-y,
第三個元素2,你把加法套到這三個元素當然是2+b-y

: 不過又回到第一個問題來,還是他在Block 裡能表現出特別的意義來?
: 三、舉的這個例子我能看懂
: 跑 g[x_]:=Block[{w=10,xx=x},2xx+w]
: g[10] 會得到30。
: 但是例子中的 {ww,xx}和 {ww,xx}={ww,xx}
: : g[10]
: : {ww,xx}
: : g[10]=30
: : {ww,xx}={ww,xx}
: 是什麼意思呢?
ww,xx都是區域變數,所以跳出g以外那根本沒有被定義
所以還是只傳回ww,xx

: 不好意思問題蠻多的 ><
: 如果能抽空幫忙回答的話,我會非常感激的
: 感謝!


--
養花種魚數月亮賞星星

http://cydye1069.blogspot.com

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

[自動轉寄] Re: [問題] mathematica 的block用法和三個問題

作者: Passions (passion)
標題: Re: [問題] mathematica 的block用法和三個問題
時間: Sat Jul 24 15:24:23 2010

※ 引述《chungyuandye (養花種魚數月亮賞星星)》之銘言:
: ※ 引述《Passions (passion)》之銘言:
: : 想請問 Block 的用法是什麼?
: : 這一段程式碼我看了好久,也有找 help 來看,但一直都看不懂
: : 可以請高手幫我解答一下嗎:
: : f[L_?NumberQ,G_?NumberQ,K_NumberQ]:=Block[{r,n,t,y},
: : {r}=NDSolve[{包含n,t,L,G的微分方程式,n[0]==K},n,{t,0,10}];
: : y=n[a] /. r;Plus @@ ((b-y)^2)]
: : ...
: : 最後目的是對 f 這個自訂函數找 FindMinimum.
: : 我的問題是:
: : 一、Block 是不是用來創造出一個環境是使得裡面的數值不受外界影響,
: : 然後指定不受影響的數為 r,n,t,y。不太能了解用在這邊的意義為何?
: : 二、思考流程是什麼?是先從 NDSolve解出數值,再丟給 Block
: : 做不知道什麼的用途,之後再傳給 f 函數?@@
: : 三、Plus @@ 這個用法是怎麼去做計算的呢,我開新頁去測試,仍然理不出頭緒
: : 請問有沒有大大能夠幫忙回答的呢?因為還蠻急的,找了不少人問都沒有結果。
: : 感謝!
: {r,n,t,y}是區域變數,只有在f這個函數才有作用,一旦跳出f它什麼都不是
: 舉個例子
: g[x_]:=Block[{w=10,xx=x},2xx+w]
: 執行一下
: g[10]
: {ww,xx}
: g[10]=30
: {ww,xx}={ww,xx}
: 不過為了怕之前L,G,K這些變數已經有定義過,建議改成這樣
: f[L_?NumberQ,G_?NumberQ,K_NumberQ]:=Block[{LL=L,GG=G,KK=K,r,n,t,y},
: {r}=NDSolve[{包含n,t,LL,GG的微分方程式,n[0]==K},n,{t,0,10}];
: y=n[a] /. r;Plus @@ ((b-y)^2)]
: @@=Apply,Plus @@ ((b-y)^2)] =>把加法套用在((b-y)^2)上

謝謝你的回信!

不過仍然有幾個問題:

一、關於 Block 的用法,仍有一點不清楚。

原式中的 Block[{A},B ; C ; D ]

分號; 在 mathematica 中,應該是代表指令結束,

但在這邊還是被包在 Block 裡面,這真的是非常奇怪,不知道語法是什麼。

二、Plus @@ ((b-y)^2)]

這個跑出來的結果會是 2 + b - y

這真的是非常奇怪,直接把它當2看,這樣不就把平方的意義給弄掉了嗎

不過又回到第一個問題來,還是他在Block 裡能表現出特別的意義來?

三、舉的這個例子我能看懂

跑 g[x_]:=Block[{w=10,xx=x},2xx+w]

g[10] 會得到30。

但是例子中的 {ww,xx}和 {ww,xx}={ww,xx}

: g[10]
: {ww,xx}
: g[10]=30
: {ww,xx}={ww,xx}

是什麼意思呢?

不好意思問題蠻多的 ><

如果能抽空幫忙回答的話,我會非常感激的

感謝!

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

2010年7月23日 星期五

[自動轉寄] (Passions) Re: [問題] mathematica 的block用法和三個問題

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [問題] mathematica 的block用法和三個問題
時間: Sat Jul 24 07:15:47 2010

※ 引述《Passions (passion)》之銘言:
: 想請問 Block 的用法是什麼?
: 這一段程式碼我看了好久,也有找 help 來看,但一直都看不懂
: 可以請高手幫我解答一下嗎:
: f[L_?NumberQ,G_?NumberQ,K_NumberQ]:=Block[{r,n,t,y},
: {r}=NDSolve[{包含n,t,L,G的微分方程式,n[0]==K},n,{t,0,10}];
: y=n[a] /. r;Plus @@ ((b-y)^2)]
: ...
: 最後目的是對 f 這個自訂函數找 FindMinimum.
: 我的問題是:
: 一、Block 是不是用來創造出一個環境是使得裡面的數值不受外界影響,
: 然後指定不受影響的數為 r,n,t,y。不太能了解用在這邊的意義為何?
: 二、思考流程是什麼?是先從 NDSolve解出數值,再丟給 Block
: 做不知道什麼的用途,之後再傳給 f 函數?@@
: 三、Plus @@ 這個用法是怎麼去做計算的呢,我開新頁去測試,仍然理不出頭緒
: 請問有沒有大大能夠幫忙回答的呢?因為還蠻急的,找了不少人問都沒有結果。
: 感謝!

{r,n,t,y}是區域變數,只有在f這個函數才有作用,一旦跳出f它什麼都不是
舉個例子
g[x_]:=Block[{w=10,xx=x},2xx+w]
執行一下

g[10]
{ww,xx}

g[10]=30
{ww,xx}={ww,xx}

不過為了怕之前L,G,K這些變數已經有定義過,建議改成這樣
f[L_?NumberQ,G_?NumberQ,K_NumberQ]:=Block[{LL=L,GG=G,KK=K,r,n,t,y},
{r}=NDSolve[{包含n,t,LL,GG的微分方程式,n[0]==K},n,{t,0,10}];
y=n[a] /. r;Plus @@ ((b-y)^2)]

@@=Apply,Plus @@ ((b-y)^2)] =>把加法套用在((b-y)^2)上


--
養花種魚數月亮賞星星

http://cydye1069.blogspot.com

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

2010年7月18日 星期日

Re: [問題] SPSS曲線迴歸繪圖

作者: chungyuandye (養花種魚數月亮看星星) 看板: Statistics
標題: Re: [問題] SPSS曲線迴歸繪圖
時間: Sun May 18 14:38:32 2008

※ 引述《pigchang (300餘敗a象棋肉腳)》之銘言:
: 我用SPSS(v.13) 的curve estimate做曲線回歸 (Quadradic)
: 已分別求出VAR1, VAR2, VAR3 對 依變數Y 的迴歸線
: 但output是三個分開的圖...
: 我想觀察三條迴歸線之間的交點
: 不知該如何將這三條迴歸線畫在同一張圖??
: (i.e. 將三張圖重疊在一起的意思)
: 謝謝指點!

SCATTER PLOT => Overlay
這樣就可以∼∼

不然直接下語法

GRAPH
/SCATTERPLOT(OVERLAY)=y WITH var1 var2 vae3 (PAIR)
/MISSING=LISTWISE.

--
養花種魚數月亮賞星星

http://chungyuandye.blogspot.com

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

Re: [問題] 殘差分析中的缺適性檢定

作者: chungyuandye (養花種魚數月亮看星星) 看板: Statistics
標題: Re: [問題] 殘差分析中的缺適性檢定
時間: Sun Apr 13 12:05:54 2008

※ 引述《annie7853 (妮妮)》之銘言:
: 缺適性檢定這部分
: 老師上課的時候在這部分講的不是很詳細
: 課本上又只有一個變數x和y
: 但是 根據學長姐考古題出現了3個變數
: 那這樣缺適性檢定是否要一個一個做呢?
: 題目如下:
: model y=Bo+B1X1+B2X2+B3X3+e
: y: 5 4 5 3 2 1
: X1: 1 0 1 0 -1 -1
: X2: 1 1 1 -1 -1 -1
: X3: 1 1 0 0 -1 -1
: 就上列資料 求Lack Of Fit Test alfa=0.05
: 麻煩大家幫我解答一下 謝謝

SPSS要計算Lack Of Fit似乎只能用GLM來計算!

把Y當依變數,X1~X3當共變數!

在選項中勾選lack of fit

Tests of Between-Subjects Effects
Dependent Variable:y
|---------------|-----------------------|--|-----------|-------|----|
|Source |Type III Sum of Squares|df|Mean Square|F |Sig.|
|---------------|-----------------------|--|-----------|-------|----|
|Corrected Model|12.800a |3 |4.267 |16.000 |.059|
|---------------|-----------------------|--|-----------|-------|----|
|Intercept |66.667 |1 |66.667 |250.000|.004|
|---------------|-----------------------|--|-----------|-------|----|
|x1 |1.800 |1 |1.800 |6.750 |.122|
|---------------|-----------------------|--|-----------|-------|----|
|x2 |.229 |1 |.229 |.857 |.452|
|---------------|-----------------------|--|-----------|-------|----|
|x3 |.050 |1 |.050 |.188 |.707|
|---------------|-----------------------|--|-----------|-------|----|
|Error |.533 |2 |.267 | | |
|---------------|-----------------------|--|-----------|-------|----|
|Total |80.000 |6 | | | |
|---------------|-----------------------|--|-----------|-------|----|
|Corrected Total|13.333 |5 | | | |
|---------------|-----------------------|--|-----------|-------|----|
a. R Squared = .960 (Adjusted R Squared = .900)

Lack of Fit Tests
Dependent Variable:y
|-----------|--------------|--|-----------|----|----|
|Source |Sum of Squares|df|Mean Square|F |Sig.|
|-----------|--------------|--|-----------|----|----|
|Lack of Fit|.033 |1 |.033 |.067|.839|
|-----------|--------------|--|-----------|----|----|
|Pure Error |.500 |1 |.500 | | |
|-----------|--------------|--|-----------|----|----|

--
養花種魚數月亮賞星星

http://chungyuandye.blogspot.com

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

Re: [問題] 為 SPSS 的 scatter plot 加上直線

作者: chungyuandye (養花種魚數月亮看星星) 看板: Statistics
標題: Re: [問題] 為 SPSS 的 scatter plot 加上直線
時間: Mon May 12 17:29:34 2008

※ 引述《magicfx (去南半球度假)》之銘言:
: SPSS 16.0
: Graph > Scatter / Dot... > Simple Scatter > Define
: 畫出圖以後
: 我想另外加直線在這張圖上
: 方程已知
: 請要如何另外加上去呢?

圖形按兩下進入編輯模式∼

滑鼠右鍵=>Add a reference line from Equation

在Custom Equation的地方自己打上方程式


--
養花種魚數月亮賞星星

http://chungyuandye.blogspot.com

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

Re: [問題] SPSS-dummy variables的設定

作者: chungyuandye (養花種魚數月亮看星星) 站內: Statistics
標題: Re: [問題] SPSS-dummy variables的設定
時間: Tue Apr 15 12:01:02 2008

※ 引述《p78 (一根稻草換千袋米)》之銘言:
: 想請問Multiple Regression SPSS 15.0 version設定的問題
: 分析統計數據方面我可以處理,但是不知道怎麼設定dummy variables
: 我對於SAS比較熟,不了解SPSS設定的方式
: gender (2個) 和 ethnicity (5個) 是dummy variables.
: 有試過用GLM,2 Way ANOVA 就不用改dummy variables
: 但是我有13 dependent variables,想一次把全部數據跑出來.

以SPSS內部範例檔Cars.sav為例

Step 1. 使用GLM中的多變量分析

Step 2. 依變數選擇mpg, engine, horse
固定因子選擇 origin, cylinder

Ste3. 在模式中選擇自訂,並將主效果移入模式。
若要考慮交互作用,則同時選擇兩個主效果在移入模式中。


Step 4. 選擇參數估計值,這樣報表中會輸出參數估計值。

圖文版
http://chungyuandye.twbbs.org/2008/04/ptt-spss-dummy-variables.html

--
養花種魚數月亮賞星星

http://chungyuandye.twbbs.org

--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 118.232.189.101
※ 編輯: chungyuandye 來自: 218.173.131.99 (06/28 07:53)

Re: [問題] [SPSS]如何畫 Chi-Square plot?

作者: chungyuandye (養花種魚數月亮看星星) 看板: Statistics
標題: Re: [問題] [SPSS]如何畫 Chi-Square plot?
時間: Sat Mar 22 13:32:52 2008

※ 引述《magicfx (去南半球度假)》之銘言:
: 找來找去
: 好像沒有這個選項
: 只有找到別人用 SAS 寫的 code
: 可是我想要用 SPSS 來畫 chi-square plot
: 請問有這類的 add-ons 可以下載嗎?

SET MXLOOP =99999.
MATRIX.
GET DATA/VARIABLES=X1 to Xp/MISSING=OMIT.
COMPUTE NR=NROW(DATA).
COMPUTE NC=NCOL(DATA).
COMPUTE MEAN=T(CSUM(DATA)/NR).
PRINT MEAN.
COMPUTE SIGMA=(SSCP(DATA)-NR*MEAN*T(MEAN))/(NR-1).
PRINT SIGMA.
COMPUTE MD=MAKE(NR,2,-999).
COMPUTE MDTEMP=MAKE(NR,1,-999).
LOOP I=1 TO NR.
- COMPUTE MD(I,1)=DATA(I,:)*INV(SIGMA)*T(DATA(I,:)).
END LOOP.
COMPUTE MD(:,2)=(RNKORDER(MD(:,1))-0.5)/NR.
SAVE MD /OUTFILE=*.
END MATRIX.
RENAME VARIABLES COL1= MD COL2= CHISQUARE.
GRAPH
/SCATTERPLOT(BIVAR)=MD WITH CHISQUARE /MISSING=LISTWISE
/TITLE= 'The chi-square plot of the ordered distances.'.

不知道有沒有錯,有錯請指教。

Reference:
http://ieeexplore.ieee.org/iel5/8766/27769/01237706.pdf?arnumber=1237706

--
養花種魚數月亮賞星星

http://www.google.com.tw/search?hl=zh-TW&q=spss+box+cox&btnG=Google


--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 118.232.189.101
※ 編輯: chungyuandye 來自: 118.232.189.101 (03/22 13:39)
※ 編輯: chungyuandye 來自: 118.232.189.101 (03/22 13:56)
推 magicfx:喔謝謝 我今晚來研究一下  03/24 19:57

Re: [問題] 複選題轉成連續變項

作者: chungyuandye (養花種魚數月亮看星星) 看板: Statistics
標題: Re: [問題] 複選題轉成連續變項
時間: Fri Mar 21 12:32:31 2008

※ 引述《cansansiow (我)》之銘言:
: 問卷題目為:在這段期間,請問您有沒有去參加下列活動?(複選)
: 1、捐錢 2、穿紅衣 3、靜坐 4、915圍城 5、天下圍攻 6、遍地開花
: 92、沒有 95、拒答
: 今欲以此為多元迴歸之依變項,參加一項得一分,兩項得兩分...以此類推,請問眾
: 強者們,要如何將此名目變項轉換成連續變項?
: 謝謝~

Q1_1=捐錢, Q1_1=1:有選;Q1_1=0:沒選
其他的依此類推

開啟一個syntax檔,

IF (Q1_95 = 0) temp=Q1_1+Q1_2+Q1_3+.....Q1_92.
EXECUTE.

工具列=>Run=>All

--
養花種魚數月亮賞星星

http://chungyuandye.blogspot.com

--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 118.232.189.101
推 cansansiow:謝謝歐  03/21 21:13

[自動轉寄] Re: [問題] [SPSS]Q-Q plot 的 X/Y 軸如何對調?

作者: chungyuandye (養花種魚數月亮看星星) 看板: Statistics
標題: Re: [問題] [SPSS]Q-Q plot 的 X/Y 軸如何對調?
時間: Sun Mar 16 17:10:52 2008

※ 引述《magicfx (去南半球度假)》之銘言:
: 我想要把 SPSS 16 產生的 Q-Q plot 的 observed value 和 expected value 對調
: SPSS 預設把 Observed value 放在 x 軸
: expected value 放在 y 軸
: 但是我想要把 Observed value 改放在 y 軸
: expected value 放在 x 軸
: 請問該怎麼做呢?

不知道其他版本可不可以,SPSS 16是可以

Step 1. 先畫出Q-Q Plot

Step 2. 在Q-Q Plot上點兩下,進入圖形編輯

Step 3. 在Y軸上點兩下,選擇滑鼠右鍵

Step 4. Transpose Chart

只是這樣作,有什麼特別意義??

--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 118.232.189.101
推 magicfx:只是要跟課本的表達方式一致罷了 , 3q  03/17 23:43

[自動轉寄] Re: [問題] 迴歸裡的多元共線

作者: chungyuandye (chungyuandye) 看板: Statistics
標題: Re: [問題] 迴歸裡的多元共線
時間: Sat Jan 28 10:58:19 2006

※ 引述《latique (Afro)》之銘言:
: ※ 引述《sevenkiller (惡狠很的一發)》之銘言:
: : 就是自變項間是否有高度相關
: : 如果自變項間高度相關的話 會影響到對迴歸係數之假設檢測
: : 實際上操作的話 SPSS所提供之collinearity的統計包括
: : Tolerance VIF和Condition Index等
: : 這些統計是有關連性的 如Tolerance與VIF就是互為倒數 如果是Tolerance越小
: : 就表示該自變項與其他自變項間之共線性越高
: : 總而言之 共線性在分析當中是低比較好 共線性一高 就代表你的迴歸有問題
: : 解釋起來就很麻煩了
: 所以當各自變項出現共線性的關係
: 表示之間有高度相關
: 對於解釋依變項也較沒有解釋力的意思嗎∼∼

Y=b0+b1*X1+b2*X2+....+bk*Xk

Xi=c0+c1*X1+c2*X2+..+c_{i-1}*X_{i-1}+c_{i+1}*X_{i+1}+..ck*Xk
這邊可以算出一個R^2

VIF_i=1/(1-R^2)

VIF_i值很大(通常以10為門檻), 表示該自變數能被其他自變數所取代,
所以它可以不需要列在解釋變數

B=Inverse(X^T*X)*X^T*Y
如果是完全共線性, 那個X^T*X這個矩陣裡至少有兩行的元素等比例
所以行列式值為0, 那麼Inverse就不存在, 因此估計出來的B也會沒意義

應該是這樣子, 有錯請指正


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

2010年7月4日 星期日

Fwd: how to tell Mathematica to display x^-1 as x^-1 and not as 1/x

---------- Forwarded message ----------
From: "Sjoerd C. de Vries" <sjoerd.c.devr...@gmail.com>
Date: 7月4日, 下午3時08分
Subject: how to tell Mathematica to display x^-1 as x^-1 and not as 1/
x
To: comp.soft-sys.math.mathematica


Hi Nasser,

I think the construction below should work. I don't have much
experience with this, but I couldn't find any side effect.

Unprotect[Power];
Power /: MakeBoxes[Power[x_, y_], StandardForm] :=
 RowBox[{MakeBoxes[x, StandardForm], "^", MakeBoxes[y, StandardForm]}]
Protect[Power];

Cheers -- Sjoerd

On Jul 3, 2:19 pm, "Nasser M. Abbasi" <n...@12000.org> wrote:

> Hello;

> When I type

> x^-1

> Mathematica replies

> 1/x

> Is there a way to tell Mathematica to *display* terms with negative
> powers as  x^-n  and not as 1/x^n ? (i.e keep the term just like it
> would appear on paper).

> Another example, suppose I type

> expr = a*x^-2 + b *x^-1 +

> The reason I want to do this, is just for display purposes. When I print
> the expression on the screen, I'd like it to look like x^-1 and not like
> 1/x (to better match how the expression look like in the textbook)

> I am trying to avoid having to convert everything to a string, and force
> the form to the way I want.  I tried Hold functions, and
> TraditionalForm, but can't get it to print as I want.

> I am hoping there is an easy trick to do this?

> thanks
> --Nasser

Fwd: Is it possible to query current plot range values (or have

---------- Forwarded message ----------
From: Bill Rowe <readn...@sbcglobal.net>
Date: 7月4日, 下午6時10分
Subject: Is it possible to query current plot range values (or have
To: comp.soft-sys.math.mathematica


On 7/4/10 at 3:10 AM, dnqu...@gmail.com (Leo Alekseyev) wrote:

>Some of my plots contain vertical lines for alignment.  To make sure
>the lines extend from the top to the bottom of the plot frame, I
>typically give those lines large values for +y and -y coordinates.
>This has an unfortunate side effect that the directive
>PlotRange->All now considers my line to be a part of the plot, and
>rescales the plot range to display it in its entirety.  Is there a
>way to (a) either make PlotRange->All ignore this line somehow or
>(b) set the +y and -y coordinates of the line to match the current
>plot range?..

Here are two ways to place a vertical line on a plot without using
large coordinate values for y:

Plot[Cosh[x], {x, -5, 5},
 Epilog -> {Red, Line@{Scaled@{.8, 0}, Scaled@{.8, 1}}}]

or

Plot[Cosh[x], {x, -5, 5}, GridLinesStyle -> Red,
 GridLines -> {{3}, None}]

(shihori) Re: [其他]mathematica計算問題

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [其他]mathematica計算問題
時間: Sun Jul 4 22:29:39 2010

※ 引述《shihori (^^)》之銘言:
: 1.題目:Find the maximum of the function
: f(x ; y) = y
: subject to the constraint 395x-100y = 20 ; x,y皆為整數
: ans: y=1169 x=296
Maximize[{y, 395 x - 100 y == 20, 200 < x < 300, y > 0,
Element[x | y, Integers]}, {x, y}]

Reduce[{395 x - 100 y == 20, 200 < x < 300, y > 0}, {x, y},
Integers] /. {x_ == z1_ -> z1, y_ == z2_ -> z2, And -> List,
Or -> List}

: 2.求 [x^2] - [x] - 6 < 0 ;[]為高斯函數
: ans: -1 ≦ x < 3

Reduce[{Floor[x^2] - Floor[x] - 6 < 0}, x, Reals]
答案應該有錯!
{Floor[x^2] - Floor[x] - 6 < 0} /. x -> 2.9

: 3.求y=ln(x) 中x的的定義域y的值域
: 想問一下程式要怎麼寫才能表示出來?
: 4.解不等式 sin(x- π/6) >cos(x) ;0≦x<2π
: 但答案是 π/3 < 4π/3 ,跟程式跑出來的答案π/3 < 2π有出入
Reduce[{Sin[x - \[Pi]/6] > Cos[x], 0 <= x < 2 Pi}, x]

--
養花種魚數月亮賞星星

http://cydye1069.blogspot.com

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