แก้สมการเชิงอนุพันธ์ (Differential Equations)ด้วย DSolve[ ]

เราสามารถใช้ DSolve[ ] แก้สมการเชิงอนุพันธ์ได้มากมาย
ยกตัวอย่างเช่น เราจะแก้สมการ ของ Simple Harmonic Motion y''[t] == -a^2y[t]
โดยเราต้องบอก DSolve[ ] ว่า ฟังก์ชันคือ y[t] และตัีวแปรคือ t
คำตอบที่ได้จะมีค่าคงที่เช่น C[1], C[2], ..., C[k] ถ้าสมการมีคำตอบหลายตัว

In[19]:=

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

Out[19]=

{{y[t] C[1] Cos[a t] + C[2] Sin[a t]}}

ถ้าเรามีค่าเริ่มต้นเราก็ใส่ไว้รวมกับสมการได้เลย

In[20]:=

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

Out[20]=

{{y[t] Cos[a t]}}

นอกจากนี้ เราสามารถใส่ Boundary Conditions ได้ด้วยเช่น y'[0]=0 และ y[1]=1

In[21]:=

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

Out[21]=

{{y[t] Cos[a t] Sec[a]}}

เราสามารถแก้หลายๆสมการพร้อมๆกันได้
ในที่นี้เราสนใจ x1[t] และ x2[t] โดยที่สมการความสัมพันธ์เป็นดังนี้
x1''[t] = -x1[t] + (x2[t]-x1[t])
x2''[t] = (x1[t]-x2[t]) - x2[t]
x1[0] = x1'[0] = 0
x2[0]=1, x2'[0]=0
เราจะเก็บคำตอบไว้ใน ตัวแปรชื่อ soln

In[22]:=

soln = DSolve[{x1''[t]  -x1[t] + (x2[t] - x1[t]), x2''[t]  (x1[t] - x2[t]) - x2[t], x1[0] 0, x1 '[0] 0, x2[0] 1, x2 '[0] 0}, {x1[t], x2[t]}, t]

Out[22]=

{{x1[t] 1/2 (Cos[t] - Cos[3^(1/2) t]), x2[t] 1/2 (Cos[t] + Cos[3^(1/2) t])}}

soln เก็บคำตอบเราไว้ดังคาด

In[23]:=

soln

Out[23]=

{{x1[t] 1/2 (Cos[t] - Cos[3^(1/2) t]), x2[t] 1/2 (Cos[t] + Cos[3^(1/2) t])}}

เราสามารถดูกราฟของคำตอบได้ด้วยคำสั่งรูปแบบ Plot[Evaluate[...]...] หรือ ParametricPlot[Evaluate[...]... ]
x1[t] /. soln หมายความว่า ถ้าพบค่า x1[t]ที่ถูกกำหนดไว้ใน soln ให้นำมาใช้  ในที่นี้คือ 1/2 (Cos[t] - Cos[3^(1/2) t])

In[24]:=

Plot[Evaluate[x1[t]/.soln], {t, 0, 50}]

[Graphics:../HTMLFiles/chapter4_56.gif]

Out[24]=

⁃Graphics⁃

In[25]:=

Plot[Evaluate[x2[t]/.soln], {t, 0, 50}]

[Graphics:../HTMLFiles/chapter4_59.gif]

Out[25]=

⁃Graphics⁃

In[26]:=

ParametricPlot[Evaluate[{x1[t], x2[t]}/.soln], {t, 0, 50}, AspectRatio  Automatic]

[Graphics:../HTMLFiles/chapter4_62.gif]

Out[26]=

⁃Graphics⁃

ตัวอย่างต่อไปเราจะแก้สมการของ projectile motion ที่ยิงจาก (0,0) ด้วยความเร็วต้น Vi และมุม theta
เราจะเก็บคำตอบไว้ใน projsoln

In[27]:=

projsoln = DSolve[{y''[t]  -g, x''[t] == 0, x[0] 0, y[0] 0, x '[0] Cos[theta] Vi, y '[0] Sin[theta] Vi}, {y[t], x[t]}, t]

Out[27]=

{{y[t] 1/2 (-g t^2 + 2 t Vi Sin[theta]), x[t] t Vi Cos[theta]}}

แทนค่า Vi = 100, theta = Pi/4, g = 9.8

In[28]:=

RowBox[{projsoln, , /., RowBox[{{, RowBox[{Vi100, ,, , theta  Pi/4, ,, RowBox[{g, , 9.8}]}], }}]}]

Out[28]=

RowBox[{{, RowBox[{{, RowBox[{RowBox[{y[t], , RowBox[{1/2, , RowBox[{(, RowBox[{100 2^(1/2) t, -, RowBox[{9.8, , t^2}]}], )}]}]}], ,, x[t] 50 2^(1/2) t}], }}], }}]

เราวาดกราฟ (x[t],y[t]) ด้วย ParametricPlot[ ]

In[29]:=

RowBox[{ParametricPlot, [, RowBox[{RowBox[{Evaluate, [, RowBox[{{x[t], y[t]}, /., RowBox[{(, R ... [{g, , 9.8}]}], }}]}], )}]}], ]}], ,, {t, 0, 15}, ,, AspectRatio  Automatic}], ]}]

[Graphics:../HTMLFiles/chapter4_69.gif]

Out[29]=

⁃Graphics⁃

Created by Mathematica  (October 4, 2005)