2012年12月11日 星期二

Re: [問題] 在XY平面上的蛇行點找所有的最低值!?

作者: chungyuandye (養花種魚數月亮賞星星) 看板: Mathematica
標題: Re: [問題] 在XY平面上的蛇行點找所有的最低值!?
時間: Tue Nov 6 00:35:04 2012

※ 引述《candy88257 (阿泰斯)》之銘言:
: 標題: [問題] 在XY平面上的蛇行點找所有的最低值!?
: 時間: Mon Nov 5 22:36:00 2012
:
: 如:
:
: 給定一連串X座標為連續的座標點:(1,1),(2,3),(3,5),(4,2),(5,3),(6,7),(7,2),(8,5),(9,6)
: 畫起來像是上下起伏的震動圖
: 請問如何很快速的在裡面找到所有的最低點?
: 如上面串列的所有最低點為: (4,2),(7,2)
: PS.
: 不要用畫圖的去判斷
: 因為我目前要判斷的點都是好幾百個一組...
:
: --
: ※ 發信站: 批踢踢實業坊(ptt.cc)
: ◆ From: 175.181.151.43
: 推 jurian0101:If用法有誤,你要多看範例文件熟悉MM的語法了。 11/05 22:44
: → candy88257:意思是可以用If寫出來!? 11/05 23:00
: → candy88257:感謝指導! 11/05 23:00
: → chungyuandye:SplitBy[SortBy[data,Last],Last] 11/05 23:04
: → chungyuandye:SplitBy[SortBy[data, Last], Last][[1]] 少打[[1]] 11/05 23:04

如果要找全域極小,又不考慮第一組資料,那下列語法就可以
SplitBy[SortBy[data[[2;;-1]], Last], Last][[1]]

: → candy88257:我想找到所有的最低點...就像上面例子最後抓到2個最低 11/05 23:31
: → candy88257:也就是所有點按X順序在XY平面連起來,找到所有的最低點 11/05 23:32
: → candy88257:其實不能講最低點,應該是相對低點 11/05 23:37

如果要找相對高低點,那Split這個函數可發揮妙用

(* 產生100組,每組200個數字的隨機資料 *)
temp = Sort[#] & /@ RandomReal[{0, 50 Pi}, {100, 200}];
data = Map[{#, Sin[#]} &, temp, {2}]

(* 畫出所有資料的序列圖 *)
Graphics[Line[#] & /@ Partition[#, 2, 1], Epilog -> {Red, Point[#]},
Frame -> True, ImageSize -> 1000] & /@ data

(* 第一組資料的圖形 *)
Graphics[Line[#] & /@ Partition[data[[1]], 2, 1],
Epilog -> {Red, Point[data[[1]]]}, Frame -> True, ImageSize -> 1000]


(* 先分析第一組,找出所有的相對低點。若該組資料的y值若上一組資料y值小, *)
(* 就歸為同一組。 所以若長度為一則表示資料在遞增,故剔除。 *)

low = Select[Split[data[[1]], #1[[2]] > #2[[2]] &],
Length[#] > 1 &][[All, -1]]

(* 同樣的方法可應用在相對高點 *)
high = Select[Split[data[[1]], #1[[2]] <= #2[[2]] &],
Length[#] > 1 &][[All, -1]]

(* 輸出所有資料的相對高低點 *)
datalow =
Select[Split[data[[#]], #1[[2]] > #2[[2]] &],
Length[#] > 1 &][[All, -1]] & /@ Range[Length@data]

datahigh =
Select[Split[data[[#]], #1[[2]] <= #2[[2]] &],
Length[#] > 1 &][[All, -1]] & /@ Range[Length@data]

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

http://chungyuandye.twbbs.org


--
※ 發信站: 批踢踢實業坊(ptt.cc)
※ 編輯: chungyuandye 來自: 218.173.133.152 (11/06 00:48)
推 candy88257:我再研究一下,感謝!!  11/07 16:36

Re: [問題] mathmatica問題

作者: chungyuandye (養花種魚數月亮賞星星) 看板: Mathematica
標題: Re: [問題] mathmatica問題
時間: Thu Nov 29 22:48:09 2012

※ 引述《abbybao (小寶)》之銘言:
: 1. http://ppt.cc/NGy7
: 想問一下我想要畫這函數但為什麼x<0的部分不會顯示呢?
: 我程式有寫錯嗎@@

n的m次方,若m為分數且n為負數時,Mathematica會設定為複數。
因為是複數,左半邊才會沒畫出來。

Power[-1, 1/3]
Power[1,1/3]

(* (-1)^(1/3),複數 *)
Power[-1, 1/3] // N

ContourPlot[1-x^(2/3)==y,{x,-5,5},{y,-5,5}]

(* 方法1,自行做調整 *)
ContourPlot[(1-y)^3==x^2,{x,-5,5},{y,-5,5}]

(* 方法2,如果你根本不想要複數,就修改Power這個函數 *)
Unprotect[Power];
Power[x_?Negative,Rational[p_,q_?OddQ]]:=(-(-x)^(1/q))^p;
Protect[Power];

Power[-1, 1/3] // N
ContourPlot[1-x^(2/3)==y,{x,-5,5},{y,-5,5}]

: 2. mathmatica有可以找反函數的功能嗎?我打關鍵字inverse function
: 好像只能
: in: InverseFunction[f][x-7]
: out: f(-1)[x-7]
: 但我希望它是能顯示真值而不是以含數形式表示
: ex: f(x)=x-7 f(-1)(x)=x+7 <----希望顯示的值
: 不知道能不能做的到
: 拜託大家了,謝謝^^

f[x_] = x - 7;
InverseFunction[f][x]


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

http://chungyuandye.twbbs.org


--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 218.173.135.172
推 abbybao:感謝回應,不過有點疑問的地方,就是我算求雙曲線sinhx的  11/30 15:35
→ abbybao:的反函數時,差了一個正負號,不知道是甚麼原因?如圖  11/30 15:36
→ abbybao:http://tinyurl.com/cjepw9e  11/30 15:37
推 jurian0101:第九版多了一個Surd[R,n]函數,求實數n次方實根  11/30 18:04
→ abbybao:可惜還沒有第九版,過一陣子論壇應該就有破解板的了  12/02 12:39

Re: [問題] mathematica可以做以下運算嗎?

作者: chungyuandye (養花種魚數月亮賞星星) 看板: Mathematica
標題: Re: [問題] mathematica可以做以下運算嗎?
時間: Fri Nov 30 23:40:31 2012

※ 引述《abbybao (小寶)》之銘言:
: 1.http://tinyurl.com/cqp98qj
: 算極座標的面積

r=4-Cos[x];
2*Integrate[r^2*Pi/(2Pi),{x,0,Pi}]

: 2.http://tinyurl.com/cx9fuut
: 算卡氏座標的弧長

y=Log[1-x^2];
temp=Refine[
Simplify[Sqrt[Dt[x]^2+Dt[y]^2]]/Dt[x],{Dt[x]>0,0<x<1/4}]
ans1=Integrate[temp,{x,0,1/4}]
ans2=TrigToExp@ans1//Simplify


: 3.http://tinyurl.com/cfrosjw
: 算極座標的弧長
: 弧長部分我只有找到arclengthfactor這個指令,但是他範例只有參數座標的
: 直角跟極座標不知道怎麼寫?

r=1+Cos[x];
f[x_]={r*Cos[x],r*Sin[x]}
D[f[x],x]
Integrate[Sqrt[D[f[x],x].D[f[x],x]],{x,0,2Pi}]


: 4.http://tinyurl.com/c2mb8pd
: 算旋轉體的體積

Limit[Sum[Pi*(4/((2/n*i)^2+4))^2*2/n,{i,1,n}],n->Infinity]

Manipulate[
Graphics3D[
Cylinder[{{2/n*(#-1),0,0},{2/n*#,0,0}},
4/((2/n*#)^2+4)]&/@Range[n]],
{n,5,200,5}
]

: 5. 1/A = 1/9955 + 1/9956 + ...... + 1/9960
: 求 A 之值 ? (ANS:約1647)
: 這題我用程式寫出來卻有兩個解 π 跟 1/π
: 不知哪邊寫錯了?
: http://ppt.cc/4P7B
: 拜託大家了,謝謝^^


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

http://chungyuandye.twbbs.org


--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 218.173.135.172
→ abbybao:感恩,我弄出來了^^  12/01 12:55

Re: [問題] 矩陣行列變換

作者: chungyuandye (養花種魚數月亮賞星星) 看板: Mathematica
標題: Re: [問題] 矩陣行列變換
時間: Sat Nov 17 17:32:54 2012

※ 引述《xnicky2000 (愛談低調的地瓜)》之銘言:
: 我現在有個矩陣
: ┌ ┐
: │ a b c d │
: │ │
: │ │
: │ e f g h │
: │ │
: │ │
: │ i j k l │
: │ │
: │ │
: │ m n o p │
: │ │
: └ ┘
: 如果我要把 m n o p
: 往上移變成
: ┌ ┐
: │ a b c d │
: │ │
: │ m n o p │
: │ │
: │ e f g h │
: │ │
: │ i j k l │
: │ │
: └ ┘
: 要怎樣做啊?? 可以用TABLE做嗎?
: 因為如果之後可能又有別行要換 可以比較好了解
: 感謝
: 版本是7.0的

data = {{a, b, c, d}, {e, f, g, h}, {i, j, k, l}, {m, m, 0, p}};

Insert[RotateLeft[data[[2 ;; -1]], -1], data[[1]], 1]

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

http://chungyuandye.twbbs.org


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

Re: [問題] mathematica可以做以下運算嗎?

作者: chungyuandye (養花種魚數月亮賞星星) 看板: Mathematica
標題: Re: [問題] mathematica可以做以下運算嗎?
時間: Wed Nov 21 13:29:49 2012

※ 引述《abbybao (小寶)》之銘言:
: 1.在座標中給定三個點形成一個三角形
: 要怎麼求在這個三角形內包含的所有座標?
: ex: (0,0) (3,0) (3,3)包含了(2,1) 一點(邊界的點不算)

test[x_,y_]:=Boole[{y>0&&-3+x<0&&x-y>0}][[1]]
Select[Flatten[Table[{x,y},{x,0,3},{y,0,3}],1],test@#==1&]


: 2.http://ppt.cc/nD~K
: 求兩切線

ContourPlot[2x^3+2y^3==9x*y,{x,0,10},{y,0,10},
Epilog->{Red,PointSize[0.1],Point[{6,6}]}]

m[x_,y_]=Dt[y,x]/.Solve[Dt[2x^3+2y^3==9x*y,x],Dt[y,x]][[1]]
ans={x,y}/.Solve[{2x^3+2y^3==9x*y,y-6==m[x,y](x-6)},{x,y}]

y-#[[2]]==m@@#*(x-#[[1]])&/@ans[[1;;2]]

ContourPlot[2x^3+2y^3==9x*y,{x,0,10},{y,0,10},
Epilog->{Line[{ans[[1]],{6,6}}],Line[{ans[[2]],{6,6}}],
Red,PointSize[0.1],Point[{6,6}]}]

: 3.http://ppt.cc/ydlH
: 求面積
: 我的意思是有沒有指令是可以把邊界條件打出來它會幫你算面積
: ex: 此題邊界條件為 y= 2x-x^2 , y= x-2

Integrate[Abs[2x-x^2-(x-2)],Flatten@{x,x/.Solve[2x-x^2==(x-2),x]}]

: 4.列出100的所有質數

Select[{#,PrimeQ[#]}&/@Range[100],#[[2]]==True&][[All,1]]

: 拜託大家了,謝謝^^

http://chungyuandye.twbbs.org/2010/07/mathematica_11.html
http://chungyuandye.twbbs.org/2010/07/mathematica_18.html
--
我打研究室走過 那獨坐電腦前的容顏如苦瓜的糾結
靈感不來 長壽的煙霧不散
研究室如小小的寂寞的城 恰如商管的電梯向晚

http://chungyuandye.twbbs.org

--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 118.233.128.85
※ 編輯: chungyuandye 來自: 118.233.128.85 (11/21 13:32)
推 abbybao:哇靠,太強了,真的有辦法做出來,不過第一題我執行後是  11/21 16:49
→ abbybao:空集合,不知道是哪邊出問題 http://ppt.cc/SKtY  11/21 16:50
→ abbybao:應該想說會有一點(2,1)出來  11/21 16:50
→ chungyuandye:改成test@@#  11/21 21:05

Re: [問題] FindRoot 找不到後的後續工作

作者: chungyuandye (養花種魚數月亮賞星星) 看板: Mathematica
標題: Re: [問題] FindRoot 找不到後的後續工作
時間: Thu Nov 29 23:20:21 2012

※ 引述《ntust661 (TOEFL_5!)》之銘言:
: 我想說FindRoot動作失敗後,想要取而代之改變動作
: 但是我要如何打出當FindRoot失敗後,我改變我的動作
: EX:
: A = 1 ;
: s = FindRoot[ x^2 + x + A , {x , 1} ] ;
: 以下如果找根失敗,我就要停止動作接下來改 A = 2 , 3, 4 ...
: 我要怎麼做這個動作@@
: 我看 F1 他說找根失敗,他依舊會找出奇怪的根,我就不知道判斷的依據在哪
: 有請版友幫忙了ˊˋ

f[x_,A_]:=x^2+x+A;

NestWhile[(temp=x/.Quiet@FindRoot[f[x,#],{x,1}];
Print@{#,Abs[f[temp,#]]};
If[Abs[f[temp,#]]<10^-5,#,#-1])&,10,UnsameQ,2]

不過有警告訊息不表示說是失敗,大部分都是精確度或是newton method的條件
如果不想要這些睛告訊息,FindRoot改成Quiet@FindRoot。

Mathematica的FindRoot這個指令預設是Newton method,所以要求跟的話
最好的方式就是先Plot,找到一個好個initial。




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

http://chungyuandye.twbbs.org


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

Re: [問題] mathematica旋轉體問題

作者: chungyuandye (養花種魚數月亮賞星星) 看板: Mathematica
標題: Re: [問題] mathematica旋轉體問題
時間: Sun Dec 2 09:47:52 2012

※ 引述《abbybao (小寶)》之銘言:
: 1. http://tinyurl.com/czwfb6v
: Q1:如果我想對 y=sin(x)對x=π/2 為旋轉軸旋轉要怎麼寫呢?
: ps: 圖裡是對z軸旋轉
: Q2:如果我想要標註點(1,0,0)的位置這有辦法做的到嗎?
: ps: 圖裡是點是我自己亂標的

(* RevolutionPlot3D 的轉軸必須是有原點出發的線,可參考 RevolutionAxis設定 *)
(* 所以要對x=π/2旋轉,必須把圖形平移,之後再做調整。 *)

temp = RevolutionPlot3D[Sin[Pi/2 - x], {x, -Pi/2, 0},
BoxRatios -> {1, 1, 0.4}]

Show[ListPlot3D[temp[[1, 1]] /. {x_, y_, z_} -> {x + Pi/2, y, z},
PlotStyle -> Opacity[0.5]],
Graphics3D[{PointSize[0.025], Red, Point[{1, 0, 0}]}]]


: 2. http://ppt.cc/LO-C
: 我想要看對題目裡的拋物面跟圓柱兩者圖形取交集所產生的圖形
: 這有辦法做的到嗎?

Plot3D[(x^2 + y^2)/4, {x, -10, 10}, {y, -10, 10},
RegionFunction -> Function[{x, y}, x^2 + y^2 <= 8 y],
Filling -> Bottom]

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

http://chungyuandye.twbbs.org


--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 218.173.134.123
→ abbybao:太屌了,真的弄得出來,尤其是第二題 這對我滿有幫助的  12/02 12:40
→ abbybao:第一題我執行程式後還是沒有看到點在圖型上顯示,如圖  12/02 12:42
→ abbybao:http://ppt.cc/SrnK 我把pointsize調到30還是沒看到  12/02 12:43
→ chungyuandye:30太大了,把整個點都包住了。  12/02 15:17

(akingo14) Re: [問題] 關於動態規劃

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [問題] 關於動態規劃
時間: Mon Mar 26 21:14:15 2012

※ 引述《akingo14 (阿肯狗)》之銘言:
: 非常謝謝 感激不盡:)
: ※ 引述《chungyuandye (養花種魚數月亮賞星星)》之銘言:
: : 把preti這個陣列的第一個元素設為0
: : opt這個資料包含兩個東西,第一個最小值,第二個最小值所對定的p。
: : p = p /. opt[[2]]
: : 就是把最小值所對應的p取出來,帶回去
: : 你用迴圈,所以只是為了避免下一要要用到時變數時會重複定義,將把取消罷了。
: : 最後再把你這個階段的解丟到暫存的陣列裡
: : (*同理*)

幫你改寫一下∼

Clear["Global`*"];
timeint=200;
inven=80;
Table[vtemp[1,z]={0,0},{z,10^4}];
Table[vtemp[z,1]={0,0},{z,10^4}];
mt=0.35;

{A1,B1}={A,B}/.
FindRoot[{B*Gamma[1/A+1]==70,
B^2*(Gamma[2/A+1]-Gamma[1/A+1]^2)==100},{A,1},{B,
1}];

tempti=Table[{0},{timeint+1}];

vtemp[ti_,in_]:=
vtemp[ti,in]=Block[{prob,revcdf,p,obj,opt,ptemp},
revcdf=E^-(1/(30+0.1ti)*p)^3.25;(*Exp[-(1/B1*p)^A1]*);
prob=mt*revcdf;
obj=-(prob*(p+vtemp[ti-1,in-1][[1]])+(1-prob)*
vtemp[ti-1,in][[1]]);
opt=Quiet@FindMinimum[obj,{p,tempti[[ti,-1]]}];
ptemp=p/.opt[[2]];
tempti[[ti]]=Flatten@{tempti[[ti]],ptemp};
{-opt[[1]],ptemp}]

Table[vtemp[i,j][[1]],{i,2,timeint+1},{j,inven+1,
Max[(inven-timeint+i),2],-1}]//TableForm

Table[vtemp[i,j][[2]],{i,2,timeint+1},{j,inven+1,
Max[(inven-timeint+i),2],-1}]//TableForm

Max[tempti[[#]]]&/@Range[timeint+1]

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

http://chungyuandye.twbbs.org


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

(ChenYM) Re: [其他] Mathematica繪圖問題

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [其他] Mathematica繪圖問題
時間: Wed Jul 25 17:03:57 2012

※ 引述《ChenYM (老宅男一個)》之銘言:
: 如果想利用Mathematica繪製表現積分求面積的圖形
: 我找到的答案是Plot中用Filling
: 或是RegionPlot
: 但是該怎麼設定x的範圍呢?
: 不知道有沒有前輩可以指點一下
: 謝謝

myplot[funs_, {a_, b_}, {c_, d_}] := Show[
Plot[funs, {x, a, b}],
Plot[funs, {x, a, b}, Filling -> {1 -> {2}},
RegionFunction -> Function[{x}, c < x < d]]]

myplot[{Sin[x], Cos[x]}, {0, Pi}, {1, Pi/2}]

Manipulate[myplot[{Sin[x], Cos[x]}, {0, Pi}, {1, z}], {z, 1, Pi}]

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

http://chungyuandye.twbbs.org


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

(zxcvbnm0928) Re: [問題] 關於平行計算

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [問題] 關於平行計算
時間: Wed Aug 15 15:49:04 2012

※ 引述《zxcvbnm0928 (益)》之銘言:
: ※ 引述《chungyuandye (養花種魚數月亮賞星星)》之銘言:
: : 同時輸入兩個指令,對一個Kernel來說,還是FIFO,照順序跑
: : Mathematica並不會啟動多個核心同時運算這兩個指令
: : 如果要用平行運算,可以查一下 Parallelize,還有相關的指令
: : 例如:
: : ParallelTable[Plot[Sin[n*x],{x,0,Pi}],{n,1,10}]
: : 如果你的顯示卡有支援CUDA,8.0後的版本也可使用CUDAMap或相關的CUDA運算
: : good luck~
: 謝謝你的回答
: 根據我測試的結果
: ParallelTable[Plot[Sin[n*x],{x,0,Pi}],{n,1,4}]跑出結果的時間是
: 單純的Plot[Sin[2*x],{x,0,Pi}]的四倍多一點
: 這樣應該就代表沒有使用到平行計算吧??
: 我的電腦好像是四核心的
: 我的想法是說
: 看能不能我控制使用某一核心跑一個Plot指令
: 同時使用另一個核心跑另一個Plot指令
: 十二萬分的感謝

In[1]:= $Version

Out[1]= "8.0 for Mac OS X x86 (64-bit) (November 6, 2010)"

In[2]:= First@
AbsoluteTiming@
ParallelTable[Plot[Sin[n*x], {x, 0, 2 Pi}], {n, 1, 500}]

Out[2]= 6.570605

In[3]:= First@
AbsoluteTiming@Table[Plot[Sin[n*x], {x, 0, 2 Pi}], {n, 1, 500}]

Out[3]= 8.394393

平行運算並不一定會比較快,除非你的問題是互相獨立的問題,而且要很大量的問題
否則,花在傳輸的時間就划不來的。舉例來說,傳回Sin[n*x],n=1,2,....,500
用ParallelTable時間是6.5s,不用ParallelTable時間是8.39s

我想你要的應該是兩個圖,分別交給兩個Kernel做
那應該要用WaitAll,先將兩個指令就到Kernel,但他不會先計算,
要等下一個Kernel有空才會計算

myplot = {ParallelSubmit[Plot[Sin[x], {x, 0, 2 Pi}]],
ParallelSubmit[Plot[Cos[x], {x, 0, 2 Pi}]]};

開始計算

WaitAll[myplot]

不過我覺得你的問題應該不需要用到平行院算,應該只是你函數定義的問題而已!
否則在Plot中加一個Evaluate,要求在Plot中將function先計算好在繪圖

Plot[Evaluate@Sin[x], {x, 0, Pi}]

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

http://chungyuandye.twbbs.org


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

(chopriabin) Re: [問題] u(x,t) vs x and t 的圖

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [問題] u(x,t) vs x and t 的圖
時間: Tue Oct 2 23:49:25 2012

※ 引述《chopriabin ()》之銘言:
: 想請問一下
: 想用Mathematica畫出 u(x,t) 對 x 和 t 的圖形
: u(x,t) 是這樣定義的:
: u(x,t)=max { [ (e^t) / f(t) ] [ 1-(cosh(0.5x)-1) / g(t) ] , 0 },
: where
: g(t) is defined by
: Integrate[ sqrt[s(s+2)], {s,1, g(t)} ]=e^t-1,
: and
: f(t)=e^t+0.5* Integrate[1/g(s), {s,0, e^t-1} ]
: 故u(x,t) 是由 g(t) and f(t) 組成
: 而 g(t) and f(t) 都是implicitly defined
: 這樣用Mathematica應該要用怎樣的指令 可以畫出 u(x,t) 對 x 和 t 的圖形呢?
: 謝謝囉^^

(* Assume 1<t<10 *)

gtemp[t_]:=Block[{s,x,temp,temp1,temp2},
temp=Integrate[Sqrt[s(s+2)],s];
temp1=temp/.s->1;
temp2=temp/.s->x;
{t,x}/.Quiet@FindRoot[temp2-temp1==Exp[t]-1,{x,1}]
]

g[t_]=Interpolation[Parallelize[gtemp[#]&/@Range[1,10,0.1]]][t];

Plot[g[t],{t,1,10}]

f[t_]:=Exp[t]+0.5*NIntegrate[1/g[s],{s,0,Exp[t]-1}]

Plot[f[t],{t,1,10}]

u[x_,t_]:=Max[(Exp[t]/f[t])(1-(Cosh[0.5x]-1)/g[t]),0]

Plot3D[u[x,t],{t,1,10},{x,0,10}]

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

http://chungyuandye.twbbs.org


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

Re: [問題] 兩個曲面交線的座標點

作者: mandyyu213 (mandy)
標題: Re: [問題] 兩個曲面交線的座標點
時間: Sun Sep 30 13:19:04 2012


真的非常感謝你!
這完全就是我要的答案
太感謝你了!


※ 引述《chungyuandye (養花種魚數月亮賞星星)》之銘言:
: ※ 引述《mandyyu213 (mandy)》之銘言:
: : 我有兩個三元非線性方程式
: : 用了ContourPlot3D畫了兩個曲面 我想要得到他們的交線座標點(把它存成List)
: : 有試了幾種方法 可以在圖上只顯示他們的交線
: : 可是當我存成List時 List的點卻還是整個曲面
: : 請問有什麼方法可以只存他們的交線嗎
: : ====我用的function像這樣====
: : p1 = ContourPlot3D[{3 \[Nu]*((s2 + 0.5 s3)^2 - (s1 + 0.5 s3)^2)
: : + 2 f*Sinh[0.5 (s1 + s2 + s3)]*((1 + \[Nu])*(s1 - s2)) == 0,
: : 2*f*Cosh[0.5 (s1 + s2 + s3)] +2/3 (s1^2 - s1*s2 + s2^2
: : - s1*s3 - s2*s3 + s3^2) - 1 - f^2 == 0},
: : {s1, -1, 1}, {s2, -1, 1}, {s3, -1, 1}, Mesh -> None,
: : BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> {{Green,
: : Tube[.01]}}}, Contours -> {0}, ContourStyle -> Opacity[0]]
: : (p1畫出兩個曲面的交線)
: : p2 = Cases[p1, GraphicsComplex[pts__] :> pts, Infinity][[1]];
: : (將p1的點存到p2)
: : ListPointPlot3D[p2]
: : (將p2點畫出來 可是畫出來的是那兩個曲面)
: : ===
: : 謝謝!!
: 文章內沒設定nu和f,所以假設nu=f=1;
: \[Nu] = 1; f = 1;
: p1 = ContourPlot3D[
: Evaluate@{3 \[Nu]*((s2 + 0.5 s3)^2 - (s1 + 0.5 s3)^2) +
: 2 f*Sinh[0.5 (s1 + s2 + s3)]*((1 + \[Nu])*(s1 - s2)) == 0,
: 2*f*Cosh[0.5 (s1 + s2 + s3)] +
: 2/3 (s1^2 - s1*s2 + s2^2 - s1*s3 - s2*s3 + s3^2) - 1 - f^2 ==
: 0.5}, {s1, -1, 1}, {s2, -1, 1}, {s3, -1, 1},
: BoundaryStyle -> {{1, 2} -> {{RGBColor[0, 1, 0], PointSize[.1]}}}]
: p2 = Cases[p1, GraphicsComplex[pts__] :> pts, Infinity][[1]];
: p3 = Cases[p1, {RGBColor[0, 1, 0], PointSize[0.1], pts__} :> pts,
: Infinity][[All, 1]]
: Show[{Graphics3D@Line@p2[[p3[[1]]]], Graphics3D@Line@p2[[p3[[2]]]]}]
: p1中只是要求不顯示,不表示mathematica不會輸出。
: 所以應該要求只輸出綠色的點就好∼

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

(Sakar) Re: [請益] 尋找數列中的 peaks

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [請益] 尋找數列中的 peaks
時間: Wed Nov 28 12:50:43 2012

※ 引述《Sakar (紗卡)》之銘言:
: 我有一組光譜圖對應的data要用Mathematica處理。
: 因為先前沒有相關經驗,所以還不確定自己這樣做行不行?
: 目前的想法是先找出明顯的peak位置,但peak有好幾個,
: 有相關指令可以一下子找出所有的peak嗎?
: 寫個小程式比較data數列的前後點,也就是尋找
: x_i > x_i-1 && x_i > x_i+1
: 滿足這種條件的點。
: 有相關指令可以直接完成嗎?
: 不然只好去寫個小程式一一讀值比較了。
: 謝謝

mypeak[data_List]:=Block[{temp,high,low},
temp=If[Length[data[[1]]]>=2,data,MapIndexed[{First@#2,#1}&,data]];
high=Select[Split[temp,#2[[2]]>#1[[2]]&],Length[#]>1&][[All,-1]];
low=Select[Split[temp,#2[[2]]<#1[[2]]&],Length[#]>1&][[All,-1]];
{high,low,
ListLinePlot[temp,Frame->True,
Epilog->{PointSize[0.01],Black,Point@low,Red,Point@high}]}
]

data = MapIndexed[{First@#2, #1} &,
Accumulate[RandomReal[{-1, 1}, 100]]];

(* 相對高點, 相對低點, 圖形 *)

mypeak[data]


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

http://chungyuandye.twbbs.org


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

(agnes12) Re: [問題]Covariance不等,做discriminant analysis

作者: chungyuandye (養花種魚數月亮賞星星)
標題: Re: [問題]Covariance不等,做discriminant analysis
時間: Tue Dec 4 06:51:15 2012

※ 引述《agnes12 (agnes)》之銘言:
: Hi~
: 謝謝您的回應:)
: 我發現我的課本用的名詞有點confusing,
: 我用的課本的discriminant analysis是指canonical discriminant analysis
: 而一般(我先前google時發現)discriminant analysis
: 在我的課本用的名詞是classification analysis.
: 我現在還是疑惑我能不能做canonical discriminant analysis
: 因為data is of unequal covariance matrice ><
: ※ 引述《chungyuandye (養花種魚數月亮賞星星)》之銘言:
: : http://tinyurl.com/c5su2np
: : Quadratic Discriminant Analysis for SPSS.

I think you can analyze it with Fisher's canonical discriminant analysis/

If you try Fisher's original iris dataset, you can find that
the Covariance Matrices are Unequal.

http://tinyurl.com/ckflxg9


Discrimination methods are very much affected by the nature of
the within-group covariance matrices. The covariance matrices
are assumed equal for allgroups in the Fisher, or linear discrimination
methodology. When the assumption of homogeneity of within-groups
covariance is grossly violated, quadratic discriminant analysis is
recommended.

http://www.statistical-solutions-software.com/BMDP-documents/BMDP-5M.pdf

Linear discriminant analysis (LDA) is one approach to dimensionality reduction
that makes use of a linear transformation matrix. The widely used Fisher's
LDA is "sub-optimal" when the sample class covariance matrices are unequal,
meaning that another linear transformation exists that produces lower loss in
discrimination power.

http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=940332

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

http://chungyuandye.twbbs.org


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

2011年12月14日 星期三

Re: [問題] ListLinePlot的PlotLegend

作者: chungyuandye (養花種魚數月亮賞星星) 看板: Mathematica
標題: Re: [問題] ListLinePlot的PlotLegend
時間: Wed Dec 14 20:20:04 2011

※ 引述《enunion (珍惜)》之銘言:
: 請問有人知道再用ListLinePlot時
: 怎麼弄出一個方塊上面標示什麼顏色代表什麼意思
: 我試過PlotLegend
: 但是沒用
: 查ListLinePlot的選項也沒有PlotLegend
: 請問還有什麼方法嗎?
: 感謝~

時間很趕的話,湊合用一下∼

mylegend[plot_Graphics,legend_List]:=
Block[{p=plot,l=legend,colortemp,color,Opacitytemp,temp},
(*擷取色彩*)
colortemp=Cases[p,Hue[a_,b_,c_]:>Hue[a,b,c],Infinity];
color=colortemp[[-Length@l;;-1]];

(*配對圖例及色彩*)
temp={color[[#]],l[[#]]}&/@Range[Length@color];

(*建立圖例*)

Labeled[p,
Grid[{Graphics[{#[[1]],Thickness[0.1],Opacity[0.75],
Line[{{0,0},{1,0}}]},
ImageSize->{24,24},AspectRatio->8/24,
ImagePadding->0],#[[2]]}&/@temp],

(*指定圖例位置*)
{{Right,Top}}]
]

mylegend[ListLinePlot[Table[Accumulate[RandomReal[{-1,1},250]],{3}],
Filling->0],{"a","b","c"}]

mylegend[ListLinePlot[
Table[Accumulate[RandomReal[{-1,1},250]],{10}]],
CharacterRange["a","z"][[1;;10]]]

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

http://chungyuandye.twbbs.org


--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 218.173.134.231
推 enunion:真的可以用耶!!好強@@ 真的很感謝您和harry901大!!  12/14 21:34
推 enunion:可是標示的顏色好像怪怪的 紅色的線他會標成紫色  12/14 21:41
→ enunion:亮綠色的線會標成暗綠色 這有辦法解決嗎?  12/14 21:42
→ chungyuandye:Opacity[0.75]->Opacity[1]  12/14 21:45
→ harry901:這個太強了  12/14 21:48
推 enunion:還是一樣耶 紅色還是會標成紫色  12/14 21:51

Re: [問題] 指定變值

作者: chungyuandye (養花種魚數月亮賞星星) 看板: Mathematica
標題: Re: [問題] 指定變值
時間: Thu Dec 1 21:15:24 2011

※ 引述《qoodeva (發呆的花枝)》之銘言:
: 依邊界條件而改變的材料變數
: 第一組
: ω=6;
: λ=10;
: x=0.41;
: 第二組
: ω=7;
: λ=128;
: x=0.2;
: 第三組
: ω=9;
: λ=4800;
: x=0.1;
: =========(主程式部分)=========================================
: eqn1=6ωa x+7λb ;
: eqn2=8λa +12ωb x ;
: 將eqn1,eqn2的a,b排成係數矩陣
: 6ω x 7λ
: BB= ( )
: 8λ 12ω x
: Eigensystem[BB][[1]] (特徵值)
: {0.08,0,004}
: Eigensystem[BB][[2]] (特徵向量)
: 特徵值對應到的特徵向量
: ����������
: ↓ ↓
: {{0.008,0.004},{{0.811,0.841},{0.647,0.881}}}
: ↑ ↑
: ��������������
: 選取Eigensystem[BB][[1]]中最小值,
: 所對應到Eigensystem[BB][[2]]之特徵向量即為a,b之值。
: a=Eigensystem[BB][[2,2,1]];
: b=Eigensystem[BB][[2,2,2]];
: 再將得到的a,b值反帶回eqn1可得一條function[x]的函數xop。
: xop=eqn1/.{a->a,b->b}
: ==============================================================================
: 現在我每次改變材料變數時都可以得到一條f[x]的函數,
: 並利用這個函數作圖,最後將這三組所得到的圖形繪製一起。
: 我想請問版上高手
: 1.我只能一次一次的改變數值去得到這三條方程式,
: 是否有什麼特殊指令可以更簡潔的計算一次得到這三組函數,
: 一開始有想到用map可是它是取矩陣中的元素去運算,
: 我的目的是像Do迴圈這樣,可是我的值是只能選取特定的值。

eqn1=6\[Omega]*a*x+7\[Lambda]*b;
eqn2=8\[Lambda]*a+12\[Omega]*bx;

test[fun1_,fun2_,w_,l_,xx_]:=
Block[{a,b,BB,BE,coa,cob,eqn1temp,\[Omega],\[Lambda],x,y,
eqn11,eqn22},
eqn1temp=fun1/.x->y;
\[Omega]=w;\[Lambda]=l;x=xx;
eqn11=fun1;eqn22=fun2;
BB={Coefficient[eqn11,#]&/@{a,b},
Coefficient[eqn22,#]&/@{a,b}};
BE=Eigensystem[BB];
{coa,cob}=Sort[Transpose@BE][[1,2]];
eqn1temp/.{a->coa,b->cob}
]

data = {{6, 10, 0.41}, {7, 128, 0.2}, {9, 4800, 0.1}};

Plot[test[eqn1, eqn2, #[[1]], #[[2]], #[[3]]] & /@ data, {y, 0, 1}]


: 2.選取Eigensystem[BB][[1]]中最小值,
: 所對應到Eigensystem[BB][[2]]之特徵向量即為a,b之值。
: 要怎麼自動選取最小特徵值所對應到的特徵向量,
: 用Min[]選結果都是選到{}內最小的值=_=

先轉置,在用Sort排序,第一個就是

: 3.另外請教一下,當我使用Do迴圈所得到的大量數據,
: 是否能用成表格表示??

只有想不到,沒有做不到

myresult={#[[1]],#[[2]],#[[3]],
test[eqn1,eqn2,#[[1]],#[[2]],#[[3]]]}&/@data

TableForm[myresult,TableHeadings->{Range[Length@myresult],{"\[Omega]",
"\[Lambda]","x","xop"}}]

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

http://chungyuandye.twbbs.org


--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 218.173.131.111
推 qoodeva:謝謝老師,實在太厲害了@@  12/01 21:19
推 qoodeva:老師請教一下,當中的eqn11是指eqn1嗎??  12/01 23:10
推 LPH66:不是的 仔細看 Block 裡有宣告 eqn11 和 eqn22 為暫時變數  12/02 00:22

Re: [問題] 如何使用平行運算或加快計算時間?

作者: chungyuandye (養花種魚數月亮賞星星) 看板: Mathematica
標題: Re: [問題] 如何使用平行運算或加快計算時間?
時間: Thu Dec 1 13:01:14 2011

※ 引述《harry901 (forcing to A cup)》之銘言:
: 小弟又要來求救了@___@~ 爬了文跟網路
: 原本以為只要調整mathematica裡面的系統設定
: 就可以輕輕鬆鬆的使用多核心或連到其他電腦上進行平行運算
: 但是...嗚嗚 似乎沒那麼容易>.<~

Mathematica的Parallel Computing ToolKit預設是會開啟本機的Kernel

要利用多台電腦平行運算,那你每台電腦都要安裝LightweightGrid這個Service

如果有裝的話,會自動將網域中的電腦自動抓出來。

想看你現在有多少Kernel可以使用,Kernels[]會傳回來可使用的Kernels

要用LightWeightGrid,除非真的有大量的運算,

否則花在多台電腦的資料傳輸時間基本上划不來

若你在購買Mathematica時有加購買Premier Service,那這個是贈送的。

若找不到直接拿License到官網下載


: 爬了很多文 大多數使用平行運算是用ParallelMap來運算...

這個不一定要用Map,只是Map似乎比較直覺。

: 我天真以為 調調系統 mathematica就會聰明的自己將資料分配給多核電腦計算...
: 弄了老半天就是只有用到單核 不管在for迴圈前面加上parallelize 還是其他方法
: 都沒辦法成功平行運算>.<~ 一趟下來要四五分鐘左右....
: 由於程式的複雜度關係 沒有使用Map指令
: 小弟的程式到目前為止已經快六百行了 還有2/3要寫 寫完大約一千多行
: 程式的架構大約如下
: 主程式(Cell 1): 有函數f1,f2,f3, 變數z,t 參數a,b,c, ...
: 要計算數值重積分,還要找根 算出來的東西會用到副程式
: 副程式(Cell 2_: 約有十二個小程式 每個小程式裡面都有一個迴圈(大小不一)
: 這邊最花時間.......
: 請問 有辦法利用平行運算減少運算時間嗎?
: 實在摸不著平行運算最基本的使用概念...唉....
: 到最後 好像parallelize不能使用在for迴圈裡面的樣子
: 舉個簡單的例子好了 這樣 板上的大大比較容易知道我的問題>_<
: ===============================================================
: PB[z_, t_] = 3*t^2 - Exp[t*z];(*令一個函數PB*)
: ans = Array[anss, {10, 2}];(*令答案矩陣ans*)
: For[i = 0, i <= 10, i++;(*使用迴圈計算*)
: ans[[i, 1]] = i;
: ans[[i, 2]] = NIntegrate[NIntegrate[
: PB[z, t] - Exp[z*t] , {z, -3, 3*i}], {t, 0, 10}
: ]
: ]
: TimeUsed[]
: ans

這邊雙重積分的使用比較不好,可以先查一下Mathematica的用法

PB[z_,t_]=3*t^2-Exp[t*z];
ft[i_]=Re@Integrate[PB[z,t]-Exp[z*t],{t,0,10},{z,-3,3*i},Assumptions->i>0];

test1=AbsoluteTiming@Map[{#,N@ft[#]}&,Range[10]]

(*將函數分給各個Kernel,基本上這個指令在8.0似乎是可以省略。*)

(*但在7.0就要很仔細設定*)

DistributeDefinitions[ft,PB];

test2=AbsoluteTiming@ParallelMap[{#,N@ft[#]}&,Range[10]]
test3=AbsoluteTiming@ParallelTable[{i,N@ft[i]},{i,10}]
ans=test2[[2]]

最後是迴圈的部份,Mathematica提供蠻多的函數
可以用help查LoopingConstructs,For Do While這些函數在Mathematica不一定
是最好的選擇。程式若很複雜,需要用到很多額外的副程式,盡量用Block Module
包起來。

有錯請多指教∼

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

http://chungyuandye.twbbs.org


--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 218.173.131.111
※ 編輯: chungyuandye 來自: 218.173.131.111 (12/01 13:30)
※ 編輯: chungyuandye 來自: 218.173.131.111 (12/01 13:31)
推 harry901:謝謝老師∼ 我會吸收吸收 研發看看@@  12/01 22:49

Re: [問題] mathematica的排序與篩選

作者: chungyuandye (養花種魚數月亮賞星星) 看板: Mathematica
標題: Re: [問題] mathematica的排序與篩選
時間: Wed Dec 14 13:40:53 2011

※ 引述《harry901 (forcing to A cup)》之銘言:
: 問題:mathematica有篩選與排序的語法嗎?
: 假設 我有一串列資料ANS[m,n] m列 n行 每一列的行表示不同屬性
: 現在我需要比較這些屬性依序篩選 然後排列大小 我們假設
: 人 國 英 數
: ANS[5,4] = [王小明 99 60 31]
: [李大同 75 51 66]
: [陳小偉 71 20 43]
: [高小茜 81 80 79]
: [馬英文 78 56 26]
: 若我要篩選出所有數學大於30分且國文60以上的同學,mathematica有辦法直接篩選嗎?
: 查了資料似乎沒有篩選的語法...我只能先令一個FinalANS存放篩選結果
: 步驟一,先計算出FinalANS的長度
: templength=0;
: For[i=1,i=Length[ANS],i++,
: If[ANS[[i,2]]>60&&ANS[[i,4]]>30,templength++,{}]
: ]
: 上面可以將templength計算出FinalANS所需的長度
: 步驟二,將篩選後結果放入FinalANS中,如
: FinalANS=Array[fnlans,{templength,4}];
: tempwrite=0;
: For[i=1,i=Length[ANS],i++,
: If[ANS[[i,2]]>60&&ANS[[i,4]]>30,tempwrite++;
: FinalANS[[i]]=ANS[[tempwrite]],{}]
: ]
: 這樣就可以完成篩選的答案(步驟一二應該可以合併在同一個loop裡面,一時想不出)
: 不過 像這樣的例子
: mathematica沒辦法直接用簡單的過濾方式求得嗎?
: 最後是排序FinalANS 若我要依照英文成績由高排至低
: 那麼 我又得寫類似上面的迴圈一個一個去判斷...可以用Sort等等方式慢慢擠出來
: 試問 mathematica有相關排序的語法嗎?

data = {{"王小明", 99, 60, 31},
{"李大同", 75, 51, 66}, {"陳小偉", 71, 20, 43},
{"高小茜", 81, 80, 79}, {"馬英文", 78, 56, 26}};

mydata=Select[data, #[[4]] >= 30 && #[[2]] >= 60 &]

要排序的話,可以使用 ?Sort 或是 ?SortBy 了解更多的設定

按總分
SortBy[mydata,Total]
按姓名
SortBy[mydata,First]
按國文
SortBy[mydata,#[[2]]&]
按英文
SoryBy[mydata,#[[3]]&]

Select, Cases這些都是可以作為篩選的函數,不過Pattern Matching也是ok

data /. {a_, b_, c_, d_} :> {a, b, c, d} /; b >= 60 && d >= 30

找出1,2,.....,100任意三個數加總等於100的所有組合

AbsoluteTiming@ReplaceList[Range[100],
{___,a_,___,b_,___,c_,___}:>{a,b,c}/;a+b+c==50]

AbsoluteTiming@Select[Subsets[Range[100],{3}],Total@#==50&]

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

http://chungyuandye.twbbs.org
※ 編輯: chungyuandye 來自: 218.173.131.148 (12/14 13:53)
推 harry901:真是謝謝∼ 又學到一招了...  12/14 13:55
推 harry901:篩選用Select就可以了 我居然笨到還寫迴圈  12/14 13:59

Re: [問題] ListLinePlot的PlotLegend

作者: chungyuandye (養花種魚數月亮賞星星) 看板: Mathematica
標題: Re: [問題] ListLinePlot的PlotLegend
時間: Wed Dec 14 22:31:42 2011

※ 引述《chungyuandye (養花種魚數月亮賞星星)》之銘言:
: 標題: Re: [問題] ListLinePlot的PlotLegend
: 時間: Wed Dec 14 20:20:04 2011
:
: ※ 引述《enunion (珍惜)》之銘言:
: : 請問有人知道再用ListLinePlot時
: : 怎麼弄出一個方塊上面標示什麼顏色代表什麼意思
: : 我試過PlotLegend
: : 但是沒用
: : 查ListLinePlot的選項也沒有PlotLegend
: : 請問還有什麼方法嗎?
: : 感謝~
:
: → enunion:亮綠色的線會標成暗綠色 這有辦法解決嗎? 12/14 21:42
: 推 enunion:還是一樣耶 紅色還是會標成紫色 12/14 21:51

Mathematica會自行指定顏色,若要自行指定的話,那要改一下地方

mylegend[plot_Graphics,legend_List]:=
Block[{p=plot,l=legend,colortemp,color,Opacitytemp,temp},
(*擷取色彩*)
(*判斷是否自行者定顏色*)
colortemp=If[Cases[p,RGBColor[a_,b_,c_]:>RGBColor[a,b,c],Infinity]=={},
Cases[p,Hue[a_,b_,c_]:>Hue[a,b,c],Infinity],
Cases[p,RGBColor[a_,b_,c_]:>RGBColor[a,b,c],Infinity]];
color=colortemp[[-Length@l;;-1]];

(*配對圖例及色彩*)
temp={color[[#]],l[[#]]}&/@Range[Length@color];

(*建立圖例*)Labeled[p,
Grid[{Graphics[{#[[1]],Thickness[0.1],Opacity[0.75],
Line[{{0,0},{1,0}}]},ImageSize->{24,24},
AspectRatio->8/24,ImagePadding->0],#[[2]]}&/@
temp],
(*指定圖例位置*){{Right,Top}}]
]
mylegend[ListLinePlot[Table[Accumulate[RandomReal[{-1,1},250]],{3}],
Filling->0,PlotStyle->{Red,Blue,Green}],
{"a","b","c"}]

你會發現在沒有Filling這個Option時圖形很OK,但是加上Filling後,圖形就很不協調

而且當數列多時,顏色也是很難選擇,那不如就交由系統指定來得方便!

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

http://chungyuandye.twbbs.org


--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 218.173.134.231
推 enunion:這樣就可以了∼ 不過不懂指令為什麼這樣寫XD  12/14 22:45
→ enunion:如果要自學的話 有推薦的書單嗎?  12/14 22:46

Re: [問題] 請教有關遞迴式

作者: chungyuandye (養花種魚數月亮賞星星) 看板: Mathematica
標題: Re: [問題] 請教有關遞迴式
時間: Wed Nov 16 00:02:08 2011

※ 引述《sphwny (時間過的真快>"<)》之銘言:
: 不好意思又來問問題了
: 目前小弟寫到定義分段函數
: 然後使用遞迴的方式把剩下的函數定義完
: 舉例如下:
: H[-1,x1_,x2_,x3_,n_] :=E^(-x1*x3)*(n!/x2^n)*Sqrt[1/(2*x2)]
: 然後
: H[0,x1_,x2_,x3_,n_]:=
: x3/(2*n)*H[-1,x1,x2,x3,n+1]-x1/(2*n)*H[-1,x1,x2,x3,n] /; x2==0.5*x3^2
: H[0,x1_,x2_,x3_,n_]:=
: n!*(0.5*x1*H[-1,x1,x2,x3,n-1]-0.5*x3*H[-1,x1,x2,x3,n]) /; x2!=0.5*x3^2 ; x1>0
: H[0,x1_,x2_,x3_,n_]:=
: 0.5*x1*H[-1,x1,x2,x3,n-1]-0.5x3*H[-1,x1,x2,x3,n] /; x2!=0.5*x3^2 ; x1<0
: H[0,x1_,x2_,x3_,n_]:=
: 0.5*x3*H[-1,x1,x2,x3,n] /; x2!=0.5*x3^2; x1==0
: Table[H[i,x1_,x2_,x3_,n_]:= (1/i)*H[i-2,x1,x2,x3,n+1]-
: (x3/i)*H[i-1,x1,x2,x3,n+1]-(x1/i)*H[i-1,x1,x2,x3,n], {i,12}];
: 最後再把這些函數做些計算處理的時候
: 出現了遞迴次數超過256次之類的訊息以及後面伴隨很多錯誤訊息
: 當中似乎出現了什麼問題
: 是我的定義語法錯誤了嗎?
: 還是有其他的方法來定義這樣子的遞迴函數?
: 麻煩各位大師替我解惑!! 謝謝!!


上面的函數可以改用Which試試看!至於遞迴次數可以利用$RecursionLimit=?? 來設定
上面的函數因為都是利用SetDelayed來定義,所以每次都必須重新計算,
因此當i大的時候會非常耗時,可以試著用動態定義來改善

H[-1,x1_,x2_,x3_,n_]=E^(-x1*x3)*(n!/x2^n)*Sqrt[1/(2*x2)]

H[0,x1_,x2_,x3_,n_]:=Which[
x2==0.5*x3^2,
x3/(2*n)*H[-1,x1,x2,x3,n+1]-x1/(2*n)*H[-1,x1,x2,x3,n],
x2!=0.5*x3^2&&x1>0,
n!*(0.5*x1*H[-1,x1,x2,x3,n-1]-0.5*x3*H[-1,x1,x2,x3,n]),
x2!=0.5*x3^2&&x1<0,
0.5*x1*H[-1,x1,x2,x3,n-1]-0.5x3*H[-1,x1,x2,x3,n],
x2!=0.5*x3^2&&x1==0,
0.5*x3*H[-1,x1,x2,x3,n]
]

H[i_,x1_,x2_,x3_,n_]:=H[i,x1,x2,x3,n]=(1/i)*H[i-2,x1,x2,x3,n+1]-(x3/i)*
H[i-1,x1,x2,x3,n+1]-(x1/i)*H[i-1,x1,x2,x3,n]

利用相同的原理,可以利用Fibonacci number的函數來做實驗

f1[0] = 0; f1[1] = 1; f2[0] = 0; f2[1] = 1;
f1[x_Integer] := f1[x - 1] + f1[x - 2]
f2[x_Integer] := f2[x] = f2[x - 1] + f2[x - 2]
AbsoluteTiming@f1[30]
AbsoluteTiming@f2[30]

f1沒利用動態定義,f2則是利用動態定義來定義函數。

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

http://chungyuandye.twbbs.org


--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 218.173.136.123
※ 編輯: chungyuandye 來自: 218.173.136.123 (11/16 00:04)
※ 編輯: chungyuandye 來自: 218.173.136.123 (11/16 00:10)