Categories
Mathematica

ตัวอย่างการใช้ Mathematica แก้ปัญหาแบบง่ายๆ

พอดีมีคนถามว่า
จะหาจำนวนนับที่เล็กที่สุด ที่มีคุณสมบัติครบทั้ง 2 ข้อต่อไปนี้ได้ยังไงครับ???
I
(1) ขึ้นต้นด้วยเลข 1
(2) เมื่อสลับตัวเลขหลักแรก (ซึ่งก็คือ 1) กับตัวเลขหลักสุดท้าย แล้วจำนวนใหม่มีค่าเป็น 3 เท่าของจำนวนเดิม
II
(1) ขึ้นต้นด้วยเลข 1
(2) เมื่อย้ายตัวเลขหลักแรก (ซึ่งก็คือ 1) ไปต่อหลังตัวเลขหลักสุดท้าย แล้วจำนวนใหม่มีค่าเป็น 3 เท่าของจำนวนเดิม

(http://www.pantip.com/cafe/wahkor/topic/X11362457/X11362457.html)

อันนี้เป็นคำสั่ง Mathematica แบบง่ายๆ ที่ลองเขียนเพื่อแก้ปัญหานี้
I
swapInt1[tmp_Integer] := FromDigits@ ReplacePart[IntegerDigits[tmp], {1 -> Last@IntegerDigits[tmp], Length@IntegerDigits[tmp] -> First@IntegerDigits[tmp]}]
Catch[If[swapInt1@# == 3*#, Throw[#]] & /@ Range[1, 1000000, 1]]
คำตอบที่ได้จะเป็น Null ทั้งหมด นั่นหมายความว่ามันไม่มีคำตอบ
II
swapInt2[tmp__Integer] := FromDigits@RotateLeft@IntegerDigits@tmp

Catch[If[swapInt2@# == 3*#, Throw[#]] & /@ Range[1, 1000000, 1]]

คำตอบที่ได้คือ 142857
————————————————

จาก http://www.pantip.com/cafe/wahkor/topic/X11435196/X11435196.html

เราเรียก 6 28 496 8128….
ว่าเป็นจำนวนสมบูรณ์
เพราะมันมีคุณสมบัติที่น่าสนใจคือ
6     มีตัวประกอบแท้คือ 1,2,3
6  = 1+2+3
28   มีตัวประกอบแท้คือ 1,2,4,7,14
28  = 1+2+4+7+14 28 = 1+2+3+4+5+6+7
496 มีตัวประกอบแท้คือ 1,2,4,8,16,31,62,124,248
496 = 1+2+4+8+16+31+62+124+248
496 = 1+2+3+…+31

โดยตัวเลขพวกนี้ มันมีนิยามว่า มันอยู่ในรูป 2n-1(2n-1) เมื่อ n เป็นจำนวนนับที่ทำให้  2n-1 เป็นจำนวนเฉพาะ ซึ่งจากนิยามนี้เองเราสามารถเขียนคำสั่ง Mathematica ง่ายๆเพื่อหาตัวเลขพวกนี้ได้ว่า

(# (# + 1))/2 & /@ Select[2^Range[100] – 1, PrimeQ]

(#(#+1))/2 นี้ได้มาจากการให้ (2n-1) = x –> 2n =x+1  ดังนั้น  2n-1= (x+1)/2

——————————————
จาก
http://www.pantip.com/cafe/wahkor/topic/X11494743/X11494743.html
http://www.pantip.com/cafe/wahkor/topic/X11495190/X11495190.html

มีเลขอยู่หกหลักที่ไม่เหมือนกันซักตัว
เลขชุดนี้สามารถถอดรากที่สองได้เป็นจำนวนเต็ม
เมื่อกลับเลขชุดนี้ ก็จะยังถอดรากที่สองได้เป็นจำนวนเต็มเหมือนเดิม
เช่น MINTLA —> ถอดรากได้เป็นจำนวนเต็ม       ALTNIM —> ถอดรากได้เป็นจำนวนเต็ม
ถามว่า มีจำนวนใดบ้างที่เข้ากฏนี้และมีกี่จำนวน
ตัวอย่าง Mathematica code สำหรับปัญหานี้
(* ตั้งแต่ 3 – 10 หลัก แบบไม่ซ้ำ *)
For[i = 3, i <= 10, i++,
Print[“==”, i “==”];
num = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
test = Select[Permutations[num, {i}], #[[1]] != 0 &];
For[j = 1, j <= Length@test, j++,
If[FractionalPart@N@Sqrt@FromDigits@test[[j]] == 0 &&
FractionalPart@N@Sqrt@FromDigits@Reverse@test[[j]] == 0,
Print[FromDigits@test[[j]]]]
] ]

==3 ==
169
961
==4 ==
1089
9801
==5 ==
12769
96721
==6 ==
==7 ==
1238769
9678321
==8 ==
==9 ==
==10==

—————————–
http://www.pantip.com/cafe/wahkor/topic/X11514071/X11514071.html

ABCDEF แทนจำนวนใดได้บ้าง?
กำหนดให้
(1) ABCDEF แทนจำนวนซึ่งมี 6 หลัก
(2) ตัวอักษรแต่ละตัว แทนเลขโดด 1 – 9 ซึ่งอาจซ้ำกันได้
(3) AB ถอดรากที่สองแล้วได้จำนวนเต็ม
(4) CD ถอดรากที่สองแล้วได้จำนวนเต็ม
(5) EF ถอดรากที่สองแล้วได้จำนวนเต็ม
(6) ABCDEF ถอดรากที่สองแล้วได้จำนวนเต็ม

ตัวอย่างคำสั่ง Mathematica สำหรับปัญหานี้ครับ
num = {1, 2, 3, 4, 5, 6, 7, 8, 9};
test = Tuples[num, 6];

For[j = 1, j <= Length@test, j++,
If[FractionalPart@N@Sqrt@FromDigits@test[[j, # ;; # + 1]] & /@ {1, 3, 5} == {0, 0, 0} &&
FractionalPart@N@Sqrt@FromDigits@test[[j]] == 0,
Print[FromDigits@test[[j]]]]
]

166464
646416

——–

จาก http://www.pantip.com/cafe/wahkor/topic/X12638167/X12638167.html

 

ทำด้วย Mathematica

NestList[Total[IntegerDigits[#]^3] &, ตัวเลขที่หารด้วย 3 ลงตัว, จำนวนรอบจากข้อ (2-4) ]

เช่น เริ่มที่เลข 3 จำนวน 10 รอบ NestList[Total[IntegerDigits[#]^3]&, 3, 10]

{3, 27, 351, 153, 153, 153, 153, 153, 153, 153, 153}

หรือ เริ่มที่เลข 123345 (ตามตัวอย่าง) 10 รอบ NestList[Total[IntegerDigits[#]^3]&, 123345, 10]

{123345, 252, 141, 66, 432, 99, 1458, 702, 351, 153, 153}

มันมาหยุดอยู่ที่ 153

 

Categories
Mathematica

เข้ารหัสไฟล์

สำหรับคนที่เขียน Mathematica package และต้องการเผยแพร่แต่ไม่ต้องการให้คนเห็น codes ที่ตัวเองเขียนขึ้นมา จะด้วยเหตุผลอะไรก็ตาม สามารถใช้วิธีนี้ได้ครับ ซึ่งก็คือการเข้ารหัสด้วยการใช้คำสั่ง Encode

เช่น ถ้าไฟล์ package ต้นฉบับชื่อ mysexy.m และต้องการสร้างไฟล์ที่เข้ารหัสใหม่ชื่อ newsexy.m สามารถพิมพ์ตามนี้ได้เลยครับ

Encode[“mysexy.m”,”newsexy.m”]

ซึ่งถ้าเปิดด้วยคำสั่ง FilePrint ก็จะเห็น codes ที่ถูกเข้ารหัสอย่างตามตัวอย่างนี้ครับ

FilePrint[“newsexy.m”]

(*1N!*)mcmj<hTJue’P+lKh]7t>X#r/N5>m^c0Q tvNP”8dgNdX!bm^)#ab5F5″L,r/b/rQA==`2L8#%,0)TEnnSa<,P_<“bKr)o.to:aE+huM_jty,B

Categories
IT Mathematica

P. falciparum in a patient during treatment with artesunate

หลังจากที่ดูของคนอื่นมานาน ก็เลยลองทำส่งบ้าง 🙂

มันเป็นโมเดลที่จำลองว่าปริมาณของเชื้อมาลาเรียในร่างกายของคนไข้จะเปลี่ยนไปอย่างไรระหว่างการรักษาด้วยยาอาทีซูเนท

http://demonstrations.wolfram.com/AModelOfPlasmodiumFalciparumPopulationDynamicsInAPatientDuri/

Categories
Mathematica

นี่แหละเสน่ห์ของการเขียนโปรแกรมแก้ปัญหา

จากคำถามที่ http://forums.wolfram.com/student-support/topics/26414

กำหนดให้ b มีค่าระหว่าง  0 < b < 1
a1, a2…,a5  แต่ละตัวมีค่า 0 หรือ  1 เท่านั้น

โดยกำหนดให้ a1, a2…,a5 มีความสัมพันธ์คือ
Total[Table[Subscript[a,n]/(2^n),{n,1,5}]]]

คนถามต้องการให้เขียนfunction f ด้วย Mathematica เพื่อหาค่า

a1, a2…,a5 โดยต้องการตามนี้

f[b_]:= ?????

f[1/2] = {1,0,0,0,0} ; best fit of
f[0.344] = {0,1,0,1,1} ; {a1,a2,a3,a4,a5}

ไอ้เราก็เขียนตาม step วิธีคิดที่เข้าใจ

f[b_] := Module[{ai, tryall, pos},
ai = Tuples[{0, 1}, 5];
tryall = #[[1]]/2. + #[[2]]/4. + #[[3]]/8. + #[[4]]/16 + #[[5]]/ 32. & /@ ai;
pos = Nearest[tryall, b];
ai[[#]] & /@ Flatten[Position[tryall, #] & /@ pos]
];

พอมาเจอของอีกคนที่โพสท์ตอบด้วยถึงกับต้องนั่งหัวเราะเลยเพราะมันสั้นมาก

f[b_ /; 0 < b < 1] := IntegerDigits[Round[31*b], 2, 5]

สุดยอดจริงๆ

นี่แหละเสน่ห์ของการเขียนโปรแกรมแก้ปัญหา อย่าไปยึดติด..คิดนอกกรอบบ้าง 🙂

 

Categories
Mathematica

ธงชาติไทย

ลองเล่นอะไรสนุก ๆ ครับ
ไปโหลดธงชาติจาก wikipedia.org (http://en.wikipedia.org/wiki/Thai_flag)

flag = Import[
http://upload.wikimedia.org/wikipedia/commons/thumb/a/a9/Flag_of_Thailand.svg/225px-Flag_of_Thailand.svg.png“]

จากนั้นก็เอามาทำ animation โดยอาศัย function Texture (http://reference.wolfram.com/mathematica/ref/Texture.html) ซึ่งเป็นความสามารถอันใหม่ใน Mathematica 8

Categories
Mathematica

การเรียก subroutine ของ Fortran จาก Mathematica ด้วย NETLink และ DLL

สมมุติว่า subroutine (test.f90) ที่มีคือ

subroutine sum_and_difference(i,j,resultArray)
 implicit none
 integer*4, intent (in) :: i
 integer*4, intent (in) :: j
 integer*4, dimension(2), intent (out) :: resultArray
 resultArray(1) = i + j
 resultArray(2) = i - j
 return
 end subroutine sum_and_difference

compile ด้วย gfortran เพื่อสร้าง DLL

gfortran -c test.f90
gfortran -s -shared -mrtd -o test.dll test.o
จากนั้นในส่วนของ Mathematica ก็ทำตามนี้ครับ
 In[1]:= Needs["NETLink`"]

In[2]:= SetDirectory@NotebookDirectory[]
 Out[2]= "D:\\FU"

In[3]:= $pathToDLL = FileNameJoin[{Directory[], "test.dll"}]
 Out[3]= "D:\\FU\\test.dll"

In[4]:= SumAndDiff = DefineDLLFunction["sum_and_difference_", $pathToDLL, "void", {"Int32*", "Int32*", "Int32[]"}]

In[5]:= results = MakeNETObject[{0, 0}, "System.Int32[]"]

In[6]:= SumAndDiff[10, 20, results]

In[7]:= NETObjectToExpression@results
 Out[7]= {30, -10}

เพิ่มเติมเกี่ยวกับ NetLink
http://www.wolfram.com/learningcenter/tutorialcollection/NETLinkUserGuide/

แนะนำสำหรับคนใช้ Windows ครับ ให้ใช้ gfortran ผ่าน MinGW (http://www.mingw.org/)

Categories
Mathematica

การแสดงข้อความใน Pivot table ของ Excel

เจ้าหน้าที่ IT คนหนึ่งถามผมว่าจะแสดงข้อความใน pivot table ของ Excel ต้องทำอย่างไร

โดยเค้าเอาตัวอย่างมาให้ดูตามรูปด้านล่างนี้

ข้อมูลจริงจะมี label ประมาณหลักแสนและ header มีประมาณหลักพัน

และผลที่ต้องการอยากให้มันคล้ายๆกับด้านล่างนี้

แน่นอนครับว่าถามผมแบบนี้ผมก็บอกว่าไม่รู้ครับว่าทำยังไงใน Excel

แต่ถ้าทำใน Mathematica ล่ะง่ายนิดเดียวครับ 🙂

อันนี้เป็น Mathematica codes ที่ผมเขียนครับ

(* แสดงข้อมูลทั้งหมดของ id ที่ต้องการ *)

datid[id_, data_] := Module[{pos, headls}, headls = data[[1]];
pos = Union[Position[headls, “label”], Position[headls, “label”]] //
Flatten;
Select[data, #[[pos[[1]]]] == id &]];

 

(* แสดง list ของ header ของข้อมูล *)

headerls[dat_] := Module[{ls},
ls = Union@dat[[2 ;; Length@dat, 2]];
ls
];

 

(* เแสดงข้อมูลเป็นแถวๆของ id ที่ต้องการ  *)

makerow[id_, dat_] := Module[{hdls, iddat, tmp, output, i, j},
hdls = headerls[dat];
iddat = datid[id, dat];
tmp = Table[0, {Length@hdls}];
For[i = 1, i <= Length@iddat, i++,
For[j = 1, j <= Length@hdls, j++,
If[iddat[[i, 2]] == hdls[[j]],
tmp = ReplacePart[tmp, j -> iddat[[i, 3]]]]]];
{id}~Join~tmp
];

(* function หลักสำหรับอ่านไฟล์ข้อมูล และแสดงผล*)

mtr[filename_] := Module[{dat, outls, hdrls, idls},
dat = Import[filename, {“Data”, 1}];
hdrls = headerls[dat];
idls = idlist[dat];

outls = makerow[#, dat] & /@ idls;
{{“label”}~Join~hdrls}~Join~outls
];

 

(* เขียนผลลัพท์เป็น Excel *)

Export[“outfile.xlsx”,mtr[“data.xlsx”],”XLSX”]

 

 

Categories
Mathematica

การประมาณหาพื้นที่

https://www.wolframcloud.com/obj/sompob/Published/estimatearea.nb

พอดีมีคนถามว่าจะหาพื้นที่แรเงาได้อย่างไรที่

http://topicstock.pantip.com/wahkor/topicstock/2011/07/X10848195/X10848195.html

ซึ่งวิธีหาพื้นที่นี้สามารถคำนวณได้จากการอาศัยความสัมพันธ์ของพื้นที่สี่เหลี่ยมและวงกลมครับ แต่ผมจะเสนออีกวิธีหนึ่ง ซึ่งเป็นวิธีที่เรียกว่า Monte Carlo ครับ จากรูปเราอาจประมาณได้ว่า
(พื้นที่ที่ต้องการ / พื้นที่สี่เหลี่ยม(หนึ่งหน่วย) ) มีค่าประมาณ (จำนวนจุดในพื้นที่ที่ต้องการ / จำนวนจุดในพื้นที่สี่เหลี่ยม)
ดังนั้น พื้นที่ที่ต้องการ มีค่าประมาณ (จำนวนจุดในพื้นที่ที่ต้องการ / จำนวนจุดในพื้นที่สี่เหลี่ยม)

จากสูตรที่ได้นี้จะเห็นว่าถ้ายิ่งมีจำนวนจุดมากก็จะทำให้ค่าประมาณมีความละเอียดแม่นยำมากขึ้น และในตัวอย่างนี้จำนวนจุดที่ผมใช้ก็เอามาจากจำนวนจุดของภาพ(pixel)ครับ

ดังนั้นการหาจุดในพื้นที่ที่ต้องการนั้นก็เพียงนับจุดที่มีสีตรงกับสีที่ถูกทาไว้

ด้านล่างนี้เป็น Mathematica codes ที่ผมใช้ในการคำนวณครับ
โหลดภาพจากpantip.com
pic= Import[“http://topicstock.pantip.com/wahkor/topicstock/2011/07/X10848195/X10848195-0.jpg“];

crop เอาส่วนขอบออกโดยภาพที่ได้จะมีขนาด 380×380
croppic=ImageCrop[pic,{380,380}];

ดังนั้นจำนวนจุดทั้งหมดในสี่เหลี่ยมเท่ากับ 380^2
coorls =Table[{i, j}, {i, 1, 380}, {j, 1, 380}] // Flatten[#, 1] &;

จำนวนจุดที่อยู่ในพื้นที่ที่ต้องการก็คือจุดที่มีสี RGB เท่ากับ 0-0-254
ninside = Select[ImageValue[croppic, #, “Byte”] & /@ coorls, # =={ 0,0,254} &] // Length;

แต่ยังมีจุดที่อยู่บนขอบอีก
nedge=Select[ImageValue[croppic, #, “Byte”] & /@ coorls, # != {0, 0, 254} && # != {0, 0, 0} &] // Length;

ผมประมาณเอาว่าจุดที่เป็นขอบนี้เป็นส่วนของพื้นที่ที่ต้องการประมาณ 20% ดังนั้น พื้นที่ที่ต้องการประมาณ
(ninside + 0.2*nedge)/Length@coorls
0.146298

ลองเปรียบเทียบกับที่คำนวณจากสูตรโดยตรงดูนะครับ 🙂

แต่ถ้าจะใช้ Mathematica คำนวณพื้นที่ที่ว่าให้เลยก็ง่ายนิดเดียวครับ เพียงแค่ใช้คำสั่ง RegionDifference

Area@RegionDifference[Disk[{0.5, 0.5}, 0.5], Disk[{1, 0}, 1]]

0.146381

https://www.wolframcloud.com/obj/sompob/Published/estimatearea.nb

Categories
Mathematica

ทดสอบ Mathematica CDF plugin

Wolfram เพิ่งเปิดตัว plugin สำหรับแสดง CDF ไฟล์บนหน้าเวบ WordPress.

http://wordpress.org/extend/plugins/wolfram-cdf-plugin/

ใครที่มีเวบที่ใช้ WordPress อยู่ก็ลองโหลดมาใช้ดูครับ

ทดสอบ Mathematica CDF plugin

การเตรียมไฟล์ CDF สำหรับแสดงบนเวบผ่าน plugin นี้ก็ทำได้ไม่ยากครับ

พอเสร็จก็save แล้วupload ไปไว้ที่เวปเราได้เลยครับ 🙂

Categories
Mathematica

รายงานผลเลือกตั้งด้วย Mathematica

พอดีผมเห็นในทีวีหลายช่องมีการรายงานผลเลือกตั้งด้วยกราฟสวยๆเยอะมากก็คิดว่าน่าจะมาลองทำให้ดูครับว่า Mathematica ก็ทำได้ 🙂

เร็วด้วยสิ ข้อมูลที่จะเอามาเขียนกราฟผมก็ไปดูดเอามาจากเวบของคณะกรรมการการเลือกตั้งครับ

Mathematica สามารถแสดงแผนภาพหรือกราฟได้หลายแบบครับ ลองดูเพิ่มเติมได้ที่ http://reference.wolfram.com/mathematica/guide/ChartingAndInformationVisualization.html

ผมเลือกแผนภาพวงกลมครับ ซึ่งผลที่ได้ก็ตามที่แสดงด้านล่างนี้ครับ

วิธีการที่ผมทำก็คืออ่านข้อมูลจากเวบกกต http://ect.thaigov.net/election2554/report_partylist.html ด้วยคำสั่ง Import ซึ่งข้อมูลที่อ่านเข้ามาจะเป็น list ของแต่ล่ะบรรทัด ซึ่งผมดูไว้แล้วว่าข้อมูลที่จะเอามาใช้จะเริ่มที่บรรทัด 14 ถึง 53 ของข้อมูลในหน้าเวบนั้น

data=Import[“http://ect.thaigov.net/election2554/report_partylist.html”,”Data”][[14;; 53]];

ถ้าเราดูห้าบันทัดแรกของข้อมูลที่อ่านเข้ามาก็จะได้ว่า
data[[1;; 5]]
{{“หมายเลข 1″, ” พรรค เพื่อไทย”, “11,655,382 คะแนน”}, {“หมายเลข 2″,
” พรรค ชาติพัฒนาเพื่อแผ่นดิน”, “353,362 คะแนน”}, {“หมายเลข 3″,
” พรรค ประชาธิปไตยใหม่”, “103,176 คะแนน”}, {“หมายเลข 4″,
” พรรค ประชากรไทย”, “28,730 คะแนน”}, {“หมายเลข 5″,
” พรรค รักประเทศไทย”, “792,989 คะแนน”}}

จะเห็นได้ว่าข้อมูลทั้งหมดที่อ่านเข้ามาจะเป็นแบบ String ดังนั้นถ้าจะเอาข้อมูลไปใช้ในการทำแผนภูมิหรือกราฟก็ต้องมีการปรับเปลี่ยนกันหน่อย ซึ่งผมก็ได้เขียนคำสั่งง่ายๆสำหรับเปลี่ยนข้อมูลไปเป็นแบบที่ต้องการ

ber[dat_]≔Module[{tmp},tmp=Characters@dat;ToExpression@StringJoin@tmp[[9;;Length@tmp]]];
เป็นคำสั่งที่กะเอาไว้สำหรับดึงตัวแลขของหมายเลขแต่ล่ะพรรคออกมา

party[dat_] := Module[{tmp},
tmp = Characters@dat;
tmp[[6 ;; Length@tmp]] // StringJoin
];
ก็เป็นคำสั่งสำหรับดึงเอาชื่อพรรคออกมา

score[dat_] := Module[{tmp}, tmp = Characters[dat];
ToExpression@StringJoin@Select[#, # != “,” &] &@
Drop[tmp, {Length@tmp – 5, Length@tmp}]
];

ส่วนอันนี้ก็เป็นคำสั่งสำหรับดึงเอาคะแนนออกมา

เวลาจะเอาข้อมูลไปใช้ ผมก็ต้องสร้างข้อมูลขึ้นหรือจัดให้อยู่ในรูปแบบที่เราต้องการ ซึ่งในตัวอย่างผมก็จัดข้อมูลเป็นคู่ของชื่อพรรคกับคะแนนที่ได้
alldat = {party@#[[2]], score@#[[3]]} & /@ data;
อันนี้ตัวอย่างของ alldat ห้าอันแรก
alldat[[1 ;; 5]]
{{” เพื่อไทย”, 11655382}, {” ชาติพัฒนาเพื่อแผ่นดิน”,
353362}, {” ประชาธิปไตยใหม่”, 103176}, {” ประชากรไทย”,
28730}, {” รักประเทศไทย”, 792989}}

แต่ก่อนจะเอาไปใช้ผมก็ให้มันมีการเรียงลำดับจากมากไปน้อยก่อนด้วยคำสั่ง Sort

dat4plot = Sort[alldat, #1[[2]] > #2[[2]] &];

ซึ่งในรูปแผนภูมิวงกลมนั้นผมให้มันแสดงข้อมูลเป็นเปอร์เซ็นต์โดยแสดงเฉพาะพรรคที่ได้คะแนนมากกว่า 0.5%

PieChart3D[
Select[100. dat4plot[[All, 2]]/N@Total@(dat4plot[[All, 2]]), # >
0.5 &], SectorOrigin -> {Automatic, 1},
ChartLegends -> {dat4plot[[All, 1]]}]