# 工程數學期末考

### 1. 微分求極值 (20分)

* **題目內容**: 一家店在一星期可以賣出每片350元的DVD200片。假設其價格與銷售量呈線性關係，若每降價10元，可以增加銷售量20片。則他要降價多少，使收益最大？
* **解題思路**:
1. 設降價金額為 $x$ 元。
2. 降價後的價格為 $350 - x$ 元。
3. 降價 $x$ 元會使銷售量增加 $(x/10) \times 20 = 2x$ 片。
4. 降價後的銷售量為 $200 + 2x$ 片。
5. 收益函數 $R(x) = (\text{價格}) \times (\text{銷售量}) = (350 - x)(200 + 2x)$。
6. 對 $R(x)$ 求導 $R'(x)$，令 $R'(x) = 0$ 求解 $x$。
7. 驗證該點是否為最大值（例如透過二階導）。
* **Python SymPy 程式碼框架**:

    from sympy import symbols, diff, solve

    # 定義變數
    x = symbols('x')

    # 建立收益函數 R(x)
    # 價格 = 350 - x
    # 銷售量 = 200 + (x / 10) * 20
    R_x = (350 - x) * (200 + 2 * x)

    print(f"收益函數 R(x) = {R_x}")

    # 對收益函數求一階導數
    R_prime_x = diff(R_x, x)
    print(f"收益函數的一階導數 R'(x) = {R_prime_x}")

    # 令一階導數等於0，求解 x
    critical_points = solve(R_prime_x, x)
    print(f"使收益最大化的降價金額 x = {critical_points[0]}")

    # （可選）驗證二階導數確認為最大值
    # R_double_prime_x = diff(R_prime_x, x)
    # print(f"收益函數的二階導數 R''(x) = {R_double_prime_x}")
    # if R_double_prime_x.subs(x, critical_points[0]) < 0:
    #     print("該點為最大值。")

### 2. 積分求面積、體積 (20分)

* **題目內容**: 計算 $x^{2}+y^{2}+z^{2}=64$ 之球與 $x=5$ 平面形成之球帽體積。
* **解題思路**:
1. 球體半徑 $R = \sqrt{64} = 8$。
2. 球帽的高度 $h = R - 5 = 8 - 5 = 3$。
3. 使用球帽體積公式 $V = \frac{1}{3} \pi h^2 (3R - h)$。
4. 或者，透過積分圓盤面積來計算：從 $x=5$ 到 $x=R$ 對圓盤面積 $\pi (R^2 - x^2)$ 進行積分。
* **Python SymPy 程式碼框架**:

    from sympy import symbols, integrate, pi, sqrt

    # 定義變數
    x = symbols('x')

    # 球體半徑
    R = 8

    # 方法一：使用球帽體積公式
    h = R - 5
    volume_formula = (1/3) * pi * h**2 * (3*R - h)
    print(f"球帽體積 (使用公式) V = {volume_formula.evalf()}")

    # 方法二：使用積分方式計算
    # 圓盤的面積函數 A(x) = pi * (R^2 - x^2)
    area_function = pi * (R**2 - x**2)

    # 從 x=5 積分到 x=R (球體最右邊的點)
    volume_integral = integrate(area_function, (x, 5, R))
    print(f"球帽體積 (透過積分) V = {volume_integral.evalf()}")

### 3. Laplace 應用 (30分)

* **題目內容**: 使用 Laplace 轉換計算二次微分方程式之解 $\ddot{y}+2\dot{y}-3~y=2~t-e^{-2t}$ , $y(0)=1$, $\dot{y}(0)=-1$。
* **解題思路**:
1. 對微分方程兩邊取拉普拉斯變換。
* $L\{\ddot{y}\} = s^2 Y(s) - s y(0) - \dot{y}(0)$
* $L\{\dot{y}\} = s Y(s) - y(0)$
* $L\{y\} = Y(s)$
* $L\{2t\} = \frac{2}{s^2}$
* $L\{-e^{-2t}\} = -\frac{1}{s+2}$
2. 代入初始條件 $y(0)=1$, $\dot{y}(0)=-1$。
3. 解出 $Y(s)$ 的表示式。
4. 對 $Y(s)$ 進行反拉普拉斯變換，得到 $y(t)$。
* **Python SymPy 程式碼框架**:

    from sympy import symbols, Function, Eq, LaplaceTransform, exp, solve, inverse_laplace_transform
    from sympy.integrals import laplace_transform
    from sympy.abc import t, s

    # 定義函數符號和變數
    y = Function('y')
    Y_s = symbols('Y_s') # 代表 Y(s)

    # 初始條件
    y0 = 1
    ydot0 = -1

    # 定義微分方程的右側函數
    f_t = 2*t - exp(-2*t)

    # 建立拉普拉斯變換後的方程
    # (s^2 * Y_s - s * y(0) - y'(0)) + 2 * (s * Y_s - y(0)) - 3 * Y_s = L{2t} + L{-e^(-2t)}
    lhs = (s**2 * Y_s - s * y0 - ydot0) + 2 * (s * Y_s - y0) - 3 * Y_s
    rhs = 2/s**2 - 1/(s+2)

    equation = Eq(lhs, rhs)
    print(f"拉普拉斯變換後的方程: {equation}")

    # 求解 Y_s
    solution_Y_s = solve(equation, Y_s)[0]
    print(f"Y(s) 的解: {solution_Y_s}")

    # 對 Y(s) 進行反拉普拉斯變換以得到 y(t)
    y_t = inverse_laplace_transform(solution_Y_s, s, t)
    print(f"y(t) 的解: {y_t}")

### 4. Fourier 級數應用 (30分)

* **題目內容**: 將方波加在電感串聯電容並聯電阻之電路上使用傅立葉級數計算輸出電壓 Vo(電阻兩端電壓)，電路如下圖。
* **電路圖**: L1 與 (C1 並聯 R1) 串聯。輸出電壓 $V_o$ 為 R1 兩端的電壓。
* **方波函數**: $V_{i}(t)=\begin{cases}0,&0\le t<T/4\\ 1,&T/4\le t<3~T/4\end{cases}$ ,其中 $T=20e-6~sec,$ $V_{i}(t+T)=V_{i}(t)$。
* **分項計算**:
* a. $L=1e-3~H$, $C=2e-7~F$, $R=80\Omega$ 時，$\zeta=?$ $\omega=?$ $V_{o}(t)=?$ 輸出波形為何？ $V_{o}(t)$ 使用級數表示。 [cite: 2]
* b. $L=5e-4~H$, $C=1e-6~F$, $R=20\Omega$ 時，$\zeta=?$ $\omega=?$ $V_{o}(t)=?$ 輸出波形為何？ $V_{o}(t)$ 使用級數表示。
* c. $L=2e-5~H$, $C=4e-5~F$, $R=3\Omega$ 時，$\zeta=?$ $\omega=?$ $V_{o}(t)=?$ 輸出波形為何？ $V_{o}(t)$ 使用級數表示。 [cite: 3]
* 所有輸出波形時間軸 $t$ 範圍 $0 \sim 0.001~sec$。

* **解題思路**:
1. **傅立葉級數展開 $V_i(t)$**:
* 計算方波 $V_i(t)$ 的傅立葉係數 ($a_0, a_n, b_n$ 或 $C_n$)。
* 週期 $T=20e-6~sec$，基波角頻率 $\omega_0 = 2\pi/T$。
2. **計算電路的頻率響應 $H(j\omega)$**:
* 電感阻抗 $Z_L = j\omega L$。
* 電容阻抗 $Z_C = \frac{1}{j\omega C}$。
* 電阻阻抗 $Z_R = R$。
* 並聯部分的阻抗 $Z_{CR} = \frac{R \cdot Z_C}{R + Z_C} = \frac{R \cdot \frac{1}{j\omega C}}{R + \frac{1}{j\omega C}} = \frac{R}{1 + j\omega RC}$。
* 總阻抗 $Z_{total} = Z_L + Z_{CR} = j\omega L + \frac{R}{1 + j\omega RC}$。
* 傳遞函數 (頻率響應) $H(j\omega) = \frac{V_o(j\omega)}{V_i(j\omega)} = \frac{Z_{CR}}{Z_{total}}$ (分壓原理)。
3. **計算系統的固有頻率 ($\omega_n$) 和阻尼比 ($\zeta$)**:
* 將 $H(s)$ (將 $j\omega$ 替換為 $s$) 轉換為標準二階系統形式 $\frac{K\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}$，然後對應係數。
* 由 $H(s) = \frac{R_1/(L_1 R_1 C_1)}{s^2 + \frac{1}{R_1 C_1}s + \frac{1}{L_1 C_1}}$ 可得：
* $\omega_n = \frac{1}{\sqrt{L_1 C_1}}$
* $\zeta = \frac{1}{2R_1}\sqrt{\frac{L_1}{C_1}}$
4. **計算輸出電壓 $V_o(t)$ 的傅立葉級數**：
* 對於傅立葉級數的每一項，其輸出幅度是輸入幅度乘以 $|H(jn\omega_0)|$，相位是輸入相位加上 $\angle H(jn\omega_0)$。
* $V_o(t) = a_0 H(j0) + \sum_{n=1}^{\infty} (|H(jn\omega_0)| a_n \cos(n\omega_0 t + \angle H(jn\omega_0)) + |H(jn\omega_0)| b_n \sin(n\omega_0 t + \angle H(jn\omega_0)))$。

* **Python SymPy 程式碼框架**:

    from sympy import symbols, integrate, sin, cos, pi, I, exp, Abs, arg, N
    import numpy as np
    import matplotlib.pyplot as plt
    from sympy.utilities.lambdify import lambdify

    # 定義符號
    t, n, T_sym, omega_sym = symbols('t n T_sym omega_sym')
    L_sym, C_sym, R_sym = symbols('L_sym C_sym R_sym')

    # 週期與基波角頻率
    T_val = 20e-6 # sec 
    omega0 = 2 * pi / T_val # rad/s

    # --- 1. 計算輸入方波 V_i(t) 的傅立葉係數 ---
    # 方波函數 V_i(t) = 0 for 0 <= t < T/4, 1 for T/4 <= t < 3T/4 

    # 計算 a0
    a0 = (1/T_val) * integrate(1, (t, T_val/4, 3*T_val/4))
    print(f"方波傅立葉係數 a0 = {a0}")

    # 計算 an
    an_expr = (2/T_val) * integrate(1 * cos(n * omega0 * t), (t, T_val/4, 3*T_val/4))
    print(f"方波傅立葉係數 an = {an_expr}")

    # 計算 bn
    bn_expr = (2/T_val) * integrate(1 * sin(n * omega0 * t), (t, T_val/4, 3*T_val/4))
    print(f"方波傅立葉係數 bn = {bn_expr}")

    # --- 2. 計算電路的頻率響應 H(j*omega) ---
    # 傳遞函數 H(j*omega) = Z_CR / (Z_L + Z_CR)
    # Z_CR = R / (1 + j*omega*R*C)
    # Z_L = j*omega*L
    H_jw_expr = (R_sym / (1 + I * omega_sym * R_sym * C_sym)) / \
                (I * omega_sym * L_sym + R_sym / (1 + I * omega_sym * R_sym * C_sym))
    # 簡化 H(j*omega)
    H_jw_expr = H_jw_expr.simplify()
    print(f"電路傳遞函數 H(j*omega) = {H_jw_expr}")

    # --- 定義計算並繪圖的函數 ---
    def calculate_and_plot_Vo(L_val, C_val, R_val, case_name, num_terms=50): # num_terms 為傅立葉級數項數
        print(f"\n--- 情況 {case_name} ---")

        # 計算固有頻率 omega_n 和阻尼比 zeta
        omega_n_val = 1 / sqrt(L_val * C_val)
        zeta_val = 1 / (2 * R_val) * sqrt(L_val / C_val)

        print(f"  固有頻率 omega_n = {omega_n_val.evalf():.4f} rad/s")
        print(f"  阻尼比 zeta = {zeta_val.evalf():.4f}")

        # 計算輸出電壓 V_o(t) 的傅立葉級數
        Vo_t_series = a0.evalf() * Abs(H_jw_expr.subs({L_sym: L_val, C_sym: C_val, R_sym: R_val, omega_sym: 0})).evalf()

        # 累加各次諧波分量
        for n_val in range(1, num_terms + 1):
            current_omega = n_val * omega0
            H_val = H_jw_expr.subs({L_sym: L_val, C_sym: C_val, R_sym: R_val, omega_sym: current_omega})

            H_mag = Abs(H_val).evalf()
            H_phase = arg(H_val).evalf()

            an_val = an_expr.subs(n, n_val).evalf()
            bn_val = bn_expr.subs(n, n_val).evalf()

            # 傅立葉級數項疊加
            Vo_t_series += H_mag * (an_val * cos(current_omega * t + H_phase) + bn_val * sin(current_omega * t + H_phase))

        print(f"  Vo(t) (級數表示, 前 {num_terms} 項): {Vo_t_series}")

        # 繪製輸出波形
        t_values = np.linspace(0, 0.001, 500) # 時間軸 t 範圍 0~0.001sec 
        Vo_func = lambdify(t, Vo_t_series, 'numpy') # 將 SymPy 表達式轉換為 NumPy 函數
        Vo_output = Vo_func(t_values)

        plt.figure(figsize=(10, 6))
        plt.plot(t_values, Vo_output)
        plt.title(f'輸出電壓 Vo(t) for 情況 {case_name} (L={L_val}, C={C_val}, R={R_val})')
        plt.xlabel('時間 t (s)')
        plt.ylabel('電壓 Vo (V)')
        plt.grid(True)
        plt.show()

    # --- a. 參數 L=1e-3 H, C=2e-7 F, R=80 Ohm --- [cite: 2]
    calculate_and_plot_Vo(1e-3, 2e-7, 80, 'a')

    # --- b. 參數 L=5e-4 H, C=1e-6 F, R=20 Ohm --- 
    calculate_and_plot_Vo(5e-4, 1e-6, 20, 'b')

    # --- c. 參數 L=2e-5 H, C=4e-5 F, R=3 Ohm --- [cite: 3]
    calculate_and_plot_Vo(2e-5, 4e-5, 3, 'c')