Introduction to Scientific Computing 科學計算導論 Homework 1 (Due: Mar. 19, 2009)-增修版

Introduction to Scientific Computing 科學計算導論

Homework 1 (Due: Mar. 19, 2009)- 增修版

1. 考慮一初始值問題 (Initial Value Problem)

d²

dt²x =−ω0²x, 0 < t < T, (1) x(0) = x₀,

x⁰(0) = v₀, Note: ω₀ 及 T 為常數. 請選擇恰當的數值 該初始值問題的正解 (exact solution) 如下

x_exac(t) = v₀

ω₀sin ω₀t + x₀cos ω₀t.

請寫出一個 Matlab 的程式以內建的 ode45 得到 (1) 的數值解 (numerical solution):

(a) 將數值解 (xh)與正解 (xexac) 對時間作圖 (畫在同一張圖同時以不同種線表示) (b) 使用 legend 功能, 標示正解與數值解

(c) 量測數值解 (xh) 的誤差 |xh(T )− xexac(T )| (d) Question: 有辦法調整 ode45的參數使誤差變小嗎?

結果:

(a) 數值解 (xh) 與正解 (xexac)對時間作圖:

(b) 誤差:

(c) 討論: 討論一下這個小題你有什麼看法. 比方說誤差很大還是很小, 數值解和正解差很多嗎? 等等 2. 考慮一非線性初始值問題 (Initial Value Problem)

d²

dt²x =−ω0²sin x, 0 < t < T, (2) x(0) = x₀,

x⁰(0) = v₀. Note: ω₀ 及 T 為常數. 請選擇與題一相同的常數 請將第一題的程式修改對本非線性初始值問題求解.

(a) 將第一題與本題的數值解畫在同一張圖上. 改變初始條件 (選取三個不同的初始條件), 觀察兩個問 題的差異.

(b) Question: 本題的解具有週期性嗎? 如何看出 結果:

(a) 彈簧伸長量 (xspring)與單擺夾角 (xpendulum) 對時間作圖:

(b) 討論:

討論一下這個小題你有什麼看法. 比方說這兩個方程解的差異; 當初始條件 y0很小時, 兩個解相差 多少, 還有 y0很大時兩個解相差多少

3. Phase Portraits

(a) 修改第一題的程式, 將x⁰(t)對x(t)作圖 (選取五個不同的初始條件且在同一張圖上畫出).

(b) 修改第二題的程式, 將x⁰(t)對x(t)作圖 (選取五個不同的初始條件且在同一張圖上畫出).

(c) 比較兩題的 Phase Portraits 結果:

(a) Mass-Spring system:

(b) Pendulum:

(c) 討論:

討論兩個 phase portraits 的不同