แนะนำการใช้ Mathematica เบื้องต้นสำหรับนักศึกษาด้านวิทยาศาสตร์

รวบรวมโดย พงศกร สายเพ็ชร์

กลับไปหน้าแรก

ตอนที่ 1: คุณใช้ Mathematica ทำอะไรได้บ้าง

1. คุณสามารถใช้ Mathematica แทนเครื่องคิดเลขใดๆก็ได้  

เพียงแต่พิมพ์ข้อความที่ต้องการคำนวณเข้าไป แล้วกดปุ่ม Shift และ Enter ด้วยกัน ผลลัพท์ก็จะแสดงขึ้นในบรรทัดต่อไป ข้อความที่คุณพิมพ์เข้าไปจะมีคำว่า In นำหน้า ขณะที่ผลลัพท์จะมีคำว่า Out นำหน้า
ถ้ากดปุ่ม Enter อันเดียว จะขึ้นบรรทัดใหม่เท่านั้น จะไม่มีการคำนวณผลลัพท์ให้ ถ้าต้องการคำนวณผลลัพท์ของข้อความใดๆ ให้เลือกข้อความนั้นๆแล้วกดปุ่ม Shift และ Enter ด้วยกัน

ถ้าเราต้องการทราบว่า 1+1 เท่ากับอะไร เราก็พิมพ์ 1+1 แล้วกดปุ่ม Shift และ Enter ด้วยกัน ก็ได้คำตอบดังข้างล่าง

In[1]:=

1 + 1

Out[1]=

2

ถ้าเราต้องการทราบว่า 2 ยกกำลัง 3 เท่ากับอะไร เราก็พิมพ์ 2^3 แล้วกดปุ่ม Shift และ Enter ด้วยกัน

In[2]:=

2^3

Out[2]=

8

ถ้าเราต้องการทราบว่า 24 คูณกับ 25 เท่ากับอะไร เราก็พิมพ์ 24*25 หรือ 24 เว้นวรรค 25 แล้วกดปุ่ม Shift และ Enter ด้วยกัน

In[3]:=

24 * 25

Out[3]=

600

In[4]:=

24 25

Out[4]=

600

ถ้าเราต้องการทราบว่า 50 แฟคทอเรียล เท่ากับอะไร เราก็พิมพ์ 50! แล้วกดปุ่ม Shift และ Enter ด้วยกัน

In[5]:=

50 !

Out[5]=

30414093201713378043612608166064768844377641568960512000000000000

ถ้าเราต้องการทราบว่า 1 หาร 7 เท่ากับอะไร เราก็พิมพ์ 1/7 แล้วกดปุ่ม Shift และ Enter ด้วยกัน

In[6]:=

1/7

Out[6]=

1/7

คำตอบยังติดอยู่ในรูปเศษส่วน เพราะ Mathematica จะพยายามแสดงผลที่ใกล้เคียงค่าที่ถูกต้องมากที่สุดเท่าที่จะเป็นได้
ซึ่งในกรณีนี้คือจำนวนตรรกยะ 1/7
ถ้าต้องการคำตอบเป็นทศนิยม เราต้องบอก Mathematica โดยพิมพ์ 1.0/7.0 หรือ 1./7.

In[7]:=

RowBox[{1., /, 7.}]

Out[7]=

0.142857

In[8]:=

RowBox[{1., /, 7.}]

Out[8]=

0.142857

เราสามารถใช้ฟังค์ชัน N[ ] เพื่อหาค่าประมาณเป็นทศนิยม (ฟังค์ชันใน Mathematica จะอยู่ในรูป f[x,y,z,... ])

In[9]:=

N[1/7]

Out[9]=

0.142857

In[10]:=

0.142857

Out[10]=

0.142857

เราสามารถใส่จำนวนตำแหน่งนัยสำคัญที่ต้องการเข้าไปในฟังค์ชัน N[ ] ด้วย

In[11]:=

N[1/7, 2]

Out[11]=

0.14

In[12]:=

N[1/7, 20]

Out[12]=

0.14285714285714285714

In[13]:=

N[Sqrt[2], 30]

Out[13]=

1.41421356237309504880168872421

In[14]:=

N[Pi, 30]

Out[14]=

3.14159265358979323846264338328

In[15]:=

N[E, 30]

Out[15]=

2.71828182845904523536028747135

In[16]:=

N[50 !, 10]

Out[16]=

3.041409320*10^64

Mathematica รู้จักฟังค์ชันต่างมากมาย รวมถึงจำนวนเชิงซ้อน

In[17]:=

Sin[3Pi/8]

Out[17]=

Sin[(3 π)/8]

In[18]:=

RowBox[{Sin, [, RowBox[{3., , Pi/8}], ]}]

Out[18]=

0.92388

In[19]:=

Log[E]

Out[19]=

1

In[20]:=

Log[2, 16]

Out[20]=

4

In[21]:=

Exp[-1]

Out[21]=

1/

Mathematica ใช้สัญญลักษณ์ I แทนจำนวน (-1)^(1/2)

In[22]:=

RowBox[{Cos, [, RowBox[{2.5, , I}], ]}]

Out[22]=

6.13229

In[23]:=

RowBox[{Exp, [, RowBox[{1, +, RowBox[{0.5, I}]}], ]}]

Out[23]=

RowBox[{RowBox[{2.38552, }], +, RowBox[{1.30321, , }]}]

มีฟังค์ชันใน mathematical physics มากมาบ

In[24]:=

RowBox[{BesselJ, [, RowBox[{0, ,, 3.}], ]}]

Out[24]=

RowBox[{-, 0.260052}]

In[25]:=

HermiteH[9, x]

Out[25]=

30240 x - 80640 x^3 + 48384 x^5 - 9216 x^7 + 512 x^9

In[26]:=

LaguerreL[5, x]

Out[26]=

1/120 (120 - 600 x + 600 x^2 - 200 x^3 + 25 x^4 - x^5)

In[27]:=

Zeta[0]

Out[27]=

-1/2

2. คุณสามารถใช้ Mathematica สร้างกราฟทั้งสองและสามมิติ  

กราฟของ sin(x)

In[28]:=

Plot[Sin[x], {x, 0, 6Pi}]

[Graphics:HTMLFiles/chapter1_57.gif]

Out[28]=

⁃Graphics⁃

กราฟของ Damped Oscillator

In[29]:=

RowBox[{Plot, [, RowBox[{RowBox[{RowBox[{Exp, [, RowBox[{RowBox[{-, 0.2}], x}], ]}], Cos[x]}], ,, {x, 0, 10}}], ]}]

[Graphics:HTMLFiles/chapter1_60.gif]

Out[29]=

⁃Graphics⁃

กราฟของ Normal Distribution

In[30]:=

Plot[Exp[-x^2/2]/(Sqrt[2 Pi]), {x, -4, 4}]

[Graphics:HTMLFiles/chapter1_63.gif]

Out[30]=

⁃Graphics⁃

กราฟ 3 มิติ  cos(x^2+y^2)

In[31]:=

RowBox[{Plot3D, [, RowBox[{Cos[x^2 + y^2], ,, RowBox[{{, RowBox[{x, ,, RowBox[{-, 2.5}], ,, 2.5}], }}], ,, RowBox[{{, RowBox[{y, ,, RowBox[{-, 2.5}], ,, 2.5}], }}]}], ]}]

[Graphics:HTMLFiles/chapter1_66.gif]

Out[31]=

⁃SurfaceGraphics⁃

กราฟความหนาแน่น  cos(x^2+y^2) เมื่อมองจากด้านบน

In[32]:=

RowBox[{DensityPlot, [, RowBox[{Cos[x^2 + y^2], ,, RowBox[{{, RowBox[{x, ,, RowBox[{-, 2.5}], ,, 2.5}], }}], ,, RowBox[{{, RowBox[{y, ,, RowBox[{-, 2.5}], ,, 2.5}], }}]}], ]}]

[Graphics:HTMLFiles/chapter1_69.gif]

Out[32]=

⁃DensityGraphics⁃

3. คุณสามารถใช้ Mathematica ช่วยทำการคำนวณต่างๆเกี่ยวกับ Calculus  

หาลิมิตต่างๆ (พิมพ์เครื่องหมาย → โดยการพิมพ์ - แล้วพิมพ์ >)

In[33]:=

Limit[Sin[x]/x, x 0]

Out[33]=

1

In[34]:=

Limit[ArcTan[x], x Infinity]

Out[34]=

π/2

หา f '(x) เมื่อ f(x) = sin(x)

In[35]:=

D[Sin[x], x]

Out[35]=

Cos[x]

หา f '(x) เมื่อ f(x) = e^axcos(bx)

In[36]:=

D[Exp[a x] Cos[b x], x]

Out[36]=

a ^(a x) Cos[b x] - b ^(a x) Sin[b x]

หา ∫f (x) x เมื่อ f(x) = sin(ax)cos(bx)

In[37]:=

Integrate[Sin[a x] Cos[b x], x]

Out[37]=

-Cos[(a - b) x]/(2 (a - b)) - Cos[(a + b) x]/(2 (a + b))

แทนค่า a ด้วย 2 และ b ด้วย 1 โดยใช "ReplaceAll" หรือ "/." โอเปอเรเตอร์
(คุณสามารถ copy/paste ข้อความ -Cos[(a - b) x]/(2 (a - b)) - Cos[(a + b) x]/(2 (a + b))จากข้างบนได้)

In[38]:=

-Cos[(a - b) x]/(2 (a - b)) - Cos[(a + b) x]/(2 (a + b)) /. {a2, b1}

Out[38]=

-Cos[x]/2 - 1/6 Cos[3 x]

จะเกิดอะไรขึ้นถ้าเราแทนค่า a,b ด้วย 1

In[39]:=

-Cos[(a - b) x]/(2 (a - b)) - Cos[(a + b) x]/(2 (a + b)) /. {a1, b1}

1 Power :: infy : Infinite expression - encountered. More… 0

Out[39]=

ComplexInfinity

เราได้ Infinity จากการหารด้วย 0
เราจะลองหาค่าใหม่โดยให้ b=a แต่ต้น

หา ∫f (x) x เมื่อ f(x) = sin(ax)cos(ax)

In[40]:=

Integrate[Sin[a x] Cos[a x], x]

Out[40]=

-Cos[a x]^2/(2 a)

คราวนี้แทนค่า a ด้วย 1 สำเร็จ

In[41]:=

-Cos[a x]^2/(2 a) /. a1

Out[41]=

-1/2 Cos[x]^2

หา ∫_0^(2π) f(x) x เมื่อ f(x) = sin(nx)

In[42]:=

Integrate[Sin[n x], {x, 0, 2Pi}]

Out[42]=

(2 Sin[n π]^2)/n

คุณสามารถหาอนุกรมประมาณฟังค์ชันของคุณ

In[43]:=

Series[Log[1 + x], {x, 0, 5}]

Out[43]=

x - x^2/2 + x^3/3 - x^4/4 + x^5/5 + O[x]^6

In[44]:=

Series[f[x], {x, a, 8}]

Out[44]=

f[a] + f^′[a] (x - a) + 1/2 f^′′[a] (x - a)^2 + 1/6 f^(3)[a] (x - a)^3 + 1/2 ... 5 + 1/720 f^(6)[a] (x - a)^6 + (f^(7)[a] (x - a)^7)/5040 + (f^(8)[a] (x - a)^8)/40320 + O[x - a]^9

อนุกรมของ Stirling:

In[45]:=

Series[Factorial[x], {x, Infinity, 3}]

Out[45]=

^(-1 - x) x^(1 + x) ( (2 π)^(1/2) 1/x^(1/2) + 1/6  π/2^(1/2) (1/x)^(3/2) + 1/144  π/2^(1/2) (1/x)^(5/2) + O[1/x]^(7/2))

อันดับและอนุกรมก็มีให้เล่น

In[46]:=

Sum[k, {k, 1, n}]

Out[46]=

1/2 n (1 + n)

In[47]:=

Sum[k^2, {k, 1, n}]

Out[47]=

1/6 n (1 + n) (1 + 2 n)

In[48]:=

Sum[x^k/k !, {k, 0, Infinity}]

Out[48]=

^x

In[49]:=

Sum[1/(n (n + 1) (n + 2)), {n, 1, N}]

Out[49]=

(N (3 + N))/(4 (1 + N) (2 + N))

4. คุณสามารถใช้ Mathematica ช่วยเรื่องพีชคณิต

In[50]:=

Expand[(a + b)^2]

Out[50]=

a^2 + 2 a b + b^2

In[51]:=

Factor[a^2 + 2 a b + b^2]

Out[51]=

(a + b)^2

In[52]:=

Simplify[a^3 x^3 + 3 a^2 b x^3 + 3 a b^2 x^3 + b^3 x^3 - 3 a^3 x^2 y - 9 a^2 b x^2 y - 9 a b^2 ... y^2 + 9 a^2 b x y^2 + 9 a b^2 x y^2 + 3 b^3 x y^2 - a^3 y^3 - 3 a^2 b y^3 - 3 a b^2 y^3 - b^3 y^3]

Out[52]=

(a + b)^3 (x - y)^3

In[53]:=

ComplexExpand[Cosh[a + b I]]

Out[53]=

Cos[b] Cosh[a] +  Sin[b] Sinh[a]

In[54]:=

TrigExpand[Sin[ArcTan[x]]]

Out[54]=

x/(1 + x^2)^(1/2)

5. คุณสามารถใช้ Mathematica ช่วยทำแก้สมการต่าง รวมถึง differential equations

Solve[ ] สามารถแก้สมการได้หลายรูปแบบ

In[55]:=

Solve[a x^2 + b x + c  0, x]

Out[55]=

{{x (-b - (b^2 - 4 a c)^(1/2))/(2 a)}, {x (-b + (b^2 - 4 a c)^(1/2))/(2 a)}}

รวมถึงสมการที่คำตอบเป็นจำนวนเชิงซ้อน

In[56]:=

Solve[x^2 + x + 1 0, x]

Out[56]=

{{x -(-1)^(1/3)}, {x (-1)^(2/3)}}

และสมการที่มีตัวแปรหลายตัว

In[57]:=

Solve[{x + 4 y + z  1, x + 2 y + 4z 1, 4x + 2y + z1}, {x, y, z}]

Out[57]=

{{x1/8, y3/16, z1/8}}

แต่บางครั้ง Solve[ ] ก็ไม่สามารถให้คำตอบในรูปที่เราเข้าใจได้ง่ายๆ

In[58]:=

Solve[x^6 + 18x^4 + x^3 + x^2 + 10, x]

Out[58]=

{{xRoot[1 + #1^2 + #1^3 + 18 #1^4 + #1^6&, 1]}, {xRoot[1 + #1^2 + #1^3 + 1 ... #1^2 + #1^3 + 18 #1^4 + #1^6&, 5]}, {xRoot[1 + #1^2 + #1^3 + 18 #1^4 + #1^6&, 6]}}

ในกรณีเช่นนี้เราสามารถใช้ NSolve (Numerical Solve) หาคำตอบเป็นตัวเลข

In[59]:=

NSolve[x^6 + 18x^4 + x^3 + x^2 + 10, x]

Out[59]=

RowBox[{{, RowBox[{RowBox[{{, RowBox[{x, , RowBox[{RowBox[{-, 0.338255}], -, RowBox[{0 ... , RowBox[{RowBox[{0.310324, }], +, RowBox[{0.360415, , }]}]}], }}]}], }}]

นอกจากนี้เรายังสามารถแก้สมการ differential ด้วย DSolve[ ]
เราจะทดลองแก้สมการของการสั่น Simpler Harmonic Motion

In[60]:=

DSolve[y''[t]  -a y[t], y[t], t]

Out[60]=

{{y[t] C[1] Cos[a^(1/2) t] + C[2] Sin[a^(1/2) t]}}

จะเห็นว่าคำตอบข้างบนติดค่าคงที่ C1, C2 ซึ่งเราสามารถหาค่าได้ถ้าเรากำหนดเงื่อนไขเบื้องต้น (initial conditions) ให้เหมาะสม

In[61]:=

DSolve[{y''[t]  -a y[t], y[0] 1, y '[0] 0}, y[t], t]

Out[61]=

{{y[t] Cos[a^(1/2) t]}}

In[62]:=

DSolve[y''[x] - 2x y[x] + 2n y[x] 0, y[x], x]

Out[62]=

{{y[x] AiryAi[(-2 n + 2 x)/2^(2/3)] C[1] + AiryBi[(-2 n + 2 x)/2^(2/3)] C[2]}}

6. คุณสามารถใช้ Mathematica ช่วยทำการคำนวณต่างๆเกี่ยวกับ Linear Algebra  

กำหนดเวคเตอร์ด้วยลิสต์

In[63]:=

v = {1, 1, 1}

Out[63]=

{1, 1, 1}

In[64]:=

u = {1, 2, 3}

Out[64]=

{1, 2, 3}

หา dot product และ cross product

In[65]:=

Dot[u, v]

Out[65]=

6

In[66]:=

u . v

Out[66]=

6

In[67]:=

Cross[u, v]

Out[67]=

{-1, 2, -1}

บวกลบเวคเตอร์

In[68]:=

u + v

Out[68]=

{2, 3, 4}

In[69]:=

u - v

Out[69]=

{0, 1, 2}

เรากำหนดแมตริกซ์ด้วย ลิสต์ซ้อนลิสต์

In[70]:=

m = {{3, 4, 4}, {1, 2, 4}, {3, 5, 3}}

Out[70]=

{{3, 4, 4}, {1, 2, 4}, {3, 5, 3}}

แสดงให้ดูเหมือนแมตริกซ์

In[71]:=

MatrixForm[m]

Out[71]//MatrixForm=

( 3 4 4 ) 1 2 4 3 5 3

หา determinant และ inverse

In[72]:=

Det[m]

Out[72]=

-10

In[73]:=

Inverse[m]

Out[73]=

{{7/5, -4/5, -4/5}, {-9/10, 3/10, 4/5}, {1/10, 3/10, -1/5}}

In[74]:=

MatrixForm[Inverse[m]]

Out[74]//MatrixForm=

( 7 4 4 ) - -- -- 5 5 5 ... 10 10 5 1 3 1 -- -- -- 10 10 5

เราสามารถทำการคำนวณกับสัญญลักษณ์ด้วย

In[75]:=

m2 = {{a, b}, {c, d}} MatrixForm[m2]

Out[75]=

{{a, b}, {c, d}}

Out[76]//MatrixForm=

( a b ) c d

In[77]:=

MatrixForm[Inverse[m2]]

Out[77]//MatrixForm=

( d b ) ---------- ----------- -b ... c a ----------- ---------- -b c + a d -b c + a d

In[78]:=

Det[m2]

Out[78]=

-b c + a d

คุณสามารถสร้าง เวคเตอร์ และ แมตริกซ์ด้วยคำสั่ง Table[ ]

In[79]:=

vec1 = Table[Sin[(2Pi/6) x], {x, 0, 5}]

Out[79]=

{0, 3^(1/2)/2, 3^(1/2)/2, 0, -3^(1/2)/2, -3^(1/2)/2}

In[80]:=

vec2 = Table[(-1)^k, {k, 0, 5}]

Out[80]=

{1, -1, 1, -1, 1, -1}

In[81]:=

vec1 . vec2

Out[81]=

0

In[82]:=

mat1 = Table[Exp[(-1)^(m + n)], {m, 0, 5}, {n, 0, 5}]

Out[82]=

{{, 1/, , 1/, , 1/}, {1/, , 1/ ... 63309;, , 1/}, {1/, , 1/, , 1/, }}

In[83]:=

MatrixForm[mat1]

Out[83]//MatrixForm=

( 1 1 1 ) ... -- --------      

In[84]:=

mat1 . vec2

Out[84]=

{-3/ + 3 , 3/ - 3 , -3/ + 3 , 3/ - 3 , -3/ + 3 , 3/ - 3 }

In[85]:=

mat1 . vec1

Out[85]=

{0, 0, 0, 0, 0, 0}

7. คุณสามารถใช้ Mathematica เขียนโปรแกรมเพื่อแก้ปัญหาที่ซับซ้อน  

คุณสามารถเขียนโปรแกรมได้หลายรูปแบบ ทั้งแบบ procedural, functional, declarative, และ object-oriented แต่เราจะดูตัวอย่างง่ายๆก่อน

เริ่มต้นด้วยการกำหนด ฟังค์ชันของเราเอง

In[86]:=

myFunc[x_] := 2x

เราสามารถใช้ฟังค์ชันของเราในการคำนวณต่างๆ

In[87]:=

D[myFunc[x], x]

Out[87]=

2

In[88]:=

Integrate[myFunc[x], x]

Out[88]=

x^2

In[89]:=

Det[myFunc[m2]]

Out[89]=

-4 b c + 4 a d

In[90]:=

Plot[myFunc[x], {x, -1, 1}]

[Graphics:HTMLFiles/chapter1_190.gif]

Out[90]=

⁃Graphics⁃

เราสามารถกำหนด ฟังค์ชันต่างๆมาเรียงต่อกัน ให้ทำงานที่ซับซ้อนมากขึ้นเรื่อยๆ
ตัวอย่างการคำนวณ f[x+y+z] แบบอ้อมโลก:

In[91]:=

myFunc2[x_, y_, z_, f_] := Module[{a}, a = 0 ; a = a + x ; a = a + y ; a = a + z ; f[a]]

In[92]:=

myFunc2[1, 1, 1, f]

Out[92]=

f[3]

In[93]:=

myFunc2[p, q, r, Sin]

Out[93]=

Sin[p + q + r]

มี For, While loop สำหรับการทำซ้ำๆ

In[94]:=

For[i = 0, i≤5, i ++, Print["Row: ", i, "\tValues: ", i^2]]

Row: 0\tValues: 0

Row: 1\tValues: 1

Row: 2\tValues: 4

Row: 3\tValues: 9

Row: 4\tValues: 16

Row: 5\tValues: 25

In[95]:=

k = 3 ; While[k>0, Print[f[k]] ; k = k - 1]

f[3]

f[2]

f[1]

ตัวอย่างการ integrate ฟังค์ชัน f แบบสี่เหลี่ยมคางหมู โดยแบ่งช่วง (xmin, xmax) เป็น ndivisions ชิ้น

In[96]:=

trapezoidIntegrate[f_, ndivisions_Integer, {xmin_, xmax_}] := Module[{alist, h = (xmax - xmin) ... nbsp; h (-f[xmin]/2 + Apply[Plus, Map[f, alist]] - f[xmax]/2) <br />    ]

In[97]:=

trapezoidIntegrate[f, 2, {0, 1}]

Out[97]=

1/2 (f[0]/2 + f[1/2] + f[1]/2)

In[98]:=

trapezoidIntegrate[Sin, 10, {0, 1}]

Out[98]=

1/10 (Sin[1/10] + Sin[1/5] + Sin[3/10] + Sin[2/5] + Sin[1/2] + Sin[3/5] + Sin[7/10] + Sin[4/5] + Sin[9/10] + Sin[1]/2)

In[99]:=

RowBox[{trapezoidIntegrate, [, RowBox[{Sin, ,, 10, ,, RowBox[{{, RowBox[{0., ,, 1.}], }}]}], ]}]

Out[99]=

0.459315

In[100]:=

RowBox[{Integrate, [, RowBox[{Sin[x], ,, RowBox[{{, RowBox[{x, ,, 0., ,, 1.}], }}]}], ]}]

Out[100]=

0.459698

แสดงกราฟของ Hermite โพลีโนเมียลถึงอันดับ maxorder โดยที่ x อยู่ระหว่าง xmin และ xmax

In[101]:=

plotHermiteH[maxorder_, xmin_, xmax_]    := For[i = 0, i≤ maxorder, Plot[Hermi ... in, xmax}, PlotLabel  "Hermite Polynomial of Order "<>ToString[i]] ; i ++]

In[102]:=

plotHermiteH[3, -1, 1]

[Graphics:HTMLFiles/chapter1_219.gif]

[Graphics:HTMLFiles/chapter1_220.gif]

[Graphics:HTMLFiles/chapter1_221.gif]

[Graphics:HTMLFiles/chapter1_222.gif]


วันนี้ขอจบแค่นี้ก่อน  หวังว่าคุณจะพบว่า Mathematica นั้นมีประโยชน์ในการศึกษาทางด้านวิทยาศาสตร์นะครับ ตอนผมเป็นนักเรียนอยู่ Mathematica 1.0 พึ่งออกมา และผมก็ใช้ช่วยในการเรียนและทำความเข้าใจในเรื่องต่างๆได้เป็นอย่างดี

ในตอนด่อๆไป ผมจะเน้นรายละเอียดในเรื่องต่างๆมากขึ้น ถ้าใครมีอะไรสนใจเป็นพิเศษ ก็ช่วยบอกไว้ใน http://mpec.sc.mahidol.ac.th/ หรือเมล์ผม ที่ [email protected] นะครับ

Created by Mathematica  (July 18, 2005)