---
# System prepended metadata

title: HW6
tags: [CSA]

---

# HW6
[toc]

## Problem 1-(1)
### Code
1. **數據載入與模型定義**
   ```matlab
   % 載入數據
   x = [0.1, 0.15, 0.25, 0.50, 0.75, 1.0, 1.5, 3.0]; % 基質濃度 (g/L)
   mu = [0.24, 0.27, 0.34, 0.35, 0.35, 0.34, 0.33, 0.22]; % 具體生長速率 (hr^-1)
   ```
   ```matlab
   % 定義Monod模型
   monod = @(params, x) (params(1) * x) ./ (params(2) + x);
   ```
   定義Monod模型的方程，描述生長速率 (mu) 與基質濃度 (x) 之間的關係。參數 `params(1)` 表示最大生長速率 (mu_max)，`params(2)` 表示半飽和常數 (K_m)。

2. **初始猜測與模型擬合**
   ```matlab
   % 初始猜測 [mu_max, K_m]
   initial_guess = [0.4, 0.5];
   ```
   設定初始參數的猜測值，其中 `mu_max` 初始猜測值為 0.4，`K_m` 初始猜測值為 0.5。

   ```matlab
   % 使用非線性最小二乘法擬合Monod模型
   params_estimated = lsqcurvefit(monod, initial_guess, x, mu);
   ```
   我們使用 `lsqcurvefit` 函數進行非線性最小二乘擬合，來估算Monod模型的參數。

2. **提取估計參數與生成擬合曲線**
   ```matlab
   % 提取估計的參數
   mu_max_est = params_estimated(1);
   Km_est = params_estimated(2);
   ```
   提取估算出的最大生長速率 (`mu_max_est`) 和半飽和常數 (`Km_est`)。

   ```matlab
   % 生成擬合曲線
   x_fit = linspace(min(x), max(x), 100);
   mu_fit = monod(params_estimated, x_fit);
   ```
   我們生成一組基質濃度數據 `x_fit`，並使用擬合出的參數計算對應的生長速率 `mu_fit`，從而得到擬合曲線。

4. **繪圖與結果顯示**
   ```matlab
   % 繪圖
   figure; 
   plot(x, mu, 'o', 'MarkerSize', 8, 'DisplayName', 'Experimental Data');
   hold on;
   plot(x_fit, mu_fit, 'LineWidth', 2, 'DisplayName', 'Fitted Monod Model');
   xlabel('Substrate Concentration (g/L)');
   ylabel('Specific Growth Rate (\mu)');
   legend;
   title(sprintf('Fitted Monod Model\n\\mu_{max} = %.2f, K_m = %.2f', mu_max_est, Km_est));
   grid on;
   hold off;
   ```
   最後，我們繪製實驗數據點與擬合曲線，並在圖表上顯示擬合出的 `mu_max` 和 `K_m` 值。

### Result
由圖可知用Monod model擬合並不精確：

![image](https://hackmd.io/_uploads/HkOMqZESA.png =600x)

## Problem 1-(2)
### Code

1. **數據載入與轉換**
   ```matlab
   % 載入數據
   x = [0.1, 0.15, 0.25, 0.50, 0.75, 1.0, 1.5, 3.0]; % 基質濃度 (g/L)
   mu = [0.24, 0.27, 0.34, 0.35, 0.35, 0.34, 0.33, 0.22]; % 具體生長速率 (hr^-1)
   ```
   ```matlab
   % 轉換數據以用於線性模型
   Y = 1 ./ mu;
   X1 = x;
   X2 = 1 ./ x;
   X3 = ones(size(x));
   ```
   與上題不同的是，我們將數據轉換為線性模型的形式。將生長速率的倒數 (1/mu) 作為目標變量 (Y)，而基質濃度 (x)、其倒數 (1/x) 和常數項 (X3) 作為解釋變量。

2. **用線性最小二乘法求解參數**
   ```matlab
   % 用線性最小二乘法求解參數
   A = [X1', X2', X3'];
   b = Y';
   params = A \ b;
   ```
   接著使用線性最小二乘法來求解模型參數。這裡的 `A` 矩陣包含了解釋變量，而 `b` 向量包含目標變量。

3. **提取估計參數**
   ```matlab
   % 提取估計參數
   K1_mu_max = params(1);
   Km_mu_max = params(2);
   mu_max_inv = params(3);

   mu_max_est = 1 / mu_max_inv;
   K1_est = K1_mu_max * mu_max_est;
   Km_est = Km_mu_max * mu_max_est;
   ```
   從線性最小二乘法的結果中提取出估計的參數。`K1_mu_max` 和 `Km_mu_max` 分別是 `K1` 與 `mu_max` 的乘積和 `Km` 與 `mu_max` 的乘積，`mu_max_inv` 是 `mu_max` 的倒數。我們通過這些關係來計算 `mu_max`、`K1` 和 `Km` 的估計值。

4. **生成擬合曲線並繪圖**
   ```matlab
   % 生成擬合曲線
   substrate_inhibition_model = @(x) (mu_max_est * x) ./ (Km_est + x + K1_est * x.^2);
   x_fit = linspace(min(x), max(x), 100);
   mu_fit = substrate_inhibition_model(x_fit);

   % 繪圖
   figure;
   plot(x, mu, 'o', 'MarkerSize', 8, 'DisplayName', 'Experimental Data');
   hold on;
   plot(x_fit, mu_fit, 'LineWidth', 2, 'DisplayName', 'Fitted Curve');
   xlabel('Substrate Concentration (g/L)');
   ylabel('Specific Growth Rate (\mu)');
   legend;
   title(sprintf('Fitted Substrate Inhibition Model\n \\mu_{max} = %.2f, K_m = %.2f, K_1 = %.2f', mu_max_est, Km_est, K1_est));
   grid on;
   hold off;
   ```
   定義基質抑制模型方程並使用估計的參數來生成擬合曲線。然後，我們將實驗數據點和擬合曲線繪製在同一圖表中，並標註出估計的參數值。
   

### Result
![image](https://hackmd.io/_uploads/HkEBdlfHA.png =600x)


## Problem 2
### Code
1. 讀取數據
   ```matlab
   % Load data
   P = [1, 5, 10, 20, 40, 60, 100, 200, 400, 760]; % 蒸汽壓（mmHg）
   T = [-36.7, -19.6, -11.5, -2.6, 7.6, 15.4, 26.1, 42.2, 60.6, 80.1] + 273.15; % 溫度（攝氏轉為開爾文）
   ```
   這段程式碼載入了蒸汽壓 \(P\) 和對應的溫度 \(T\) 數據，並將攝氏溫度轉換為開爾文。

2. Clapeyron 方程：線性最小二乘法
   ```matlab
   %% Clapeyron Equation: log(P) = A + B/T
   % Linear Least Squares
   X1 = [ones(length(T), 1), 1./T'];
   Y1 = log10(P)';

   % Perform linear regression
   beta1 = X1 \ Y1; % Solving the normal equations
   A1 = beta1(1);
   B1 = beta1(2);

   % Model prediction
   y1_mod = 10.^(X1 * beta1);

   % Calculate the objective function (sum of squared errors)
   J1 = sum((P - y1_mod).^2);

   % Display results
   fprintf('Linear Least Squares for Clapeyron Equation:\n');
   fprintf('A = %.4f\n', A1);
   fprintf('B = %.4f\n', B1);
   fprintf('Sum of Squared Errors (J1) = %.4f\n\n', J1);
   ```
   這部分程式碼使用線性最小二乘法來擬合 Clapeyron 方程：
   1. 設置設計矩陣 \(X1\) 和響應向量 \(Y1\)。
   2. 使用矩陣求解方法得到參數 \(\beta1\)，其中包含 \(A\) 和 \(B\)。
   3. 使用模型預測值 \(y1\_mod\)。
   4. 計算目標函數（平方誤差和）。
   5. 顯示擬合結果和誤差。

3. Clapeyron 方程：非線性最小二乘法
   ```matlab
   %% Clapeyron Equation: Nonlinear Least Squares
   % Nonlinear Least Squares
   clapeyron_fun = @(params, T) 10.^(params(1) + params(2)./T);
   p0 = [0, 0]; % Initial guess for A and B
   [p, resnorm, residual] = lsqcurvefit(clapeyron_fun, p0, T, log10(P));

   % Extract parameters
   A2 = p(1);
   B2 = p(2);
   J2 = resnorm;

   % Display results
   fprintf('Nonlinear Least Squares for Clapeyron Equation:\n');
   fprintf('A = %.4f\n', A2);
   fprintf('B = %.4f\n', B2);
   fprintf('Sum of Squared Errors (J2) = %.4f\n\n', J2);
   ```
   這部分程式碼使用非線性最小二乘法來擬合 Clapeyron 方程：
   1. 定義非線性模型 `clapeyron_fun`。
   2. 設置初始參數猜測 \(p0\)。
   3. 使用 `lsqcurvefit` 函數進行擬合。
   4. 提取擬合後的參數 \(A2\) 和 \(B2\)。
   5. 計算並顯示誤差。

4. Riedel 方程：線性最小二乘法
   ```matlab
   %% Riedel Equation: log(P) = A + B/T + C*log(T) + D*T^2
   % Linear Least Squares
   X3 = [ones(length(T), 1), 1./T', log10(T)', T'.^2];
   Y3 = log10(P)';

   % Perform linear regression
   beta3 = X3 \ Y3; % Solving the normal equations
   A3 = beta3(1);
   B3 = beta3(2);
   C3 = beta3(3);
   D3 = beta3(4);

   % Model prediction
   y3_mod = 10.^(X3 * beta3);

   % Calculate the objective function (sum of squared errors)
   J3 = sum((P - y3_mod).^2);

   % Display results
   fprintf('Linear Least Squares for Riedel Equation:\n');
   fprintf('A = %.4f\n', A3);
   fprintf('B = %.4f\n', B3);
   fprintf('C = %.4f\n', C3);
   fprintf('D = %.4f\n', D3);
   fprintf('Sum of Squared Errors (J3) = %.4f\n\n', J3);
   ```
   這部分程式碼使用線性最小二乘法來擬合 Riedel 方程：
   1. 設置設計矩陣 \(X3\) 和響應向量 \(Y3\)。
   2. 使用矩陣求解方法得到參數 \(\beta3\)，其中包含 \(A\)、\(B\)、\(C\) 和 \(D\)。
   3. 使用模型預測值 \(y3\_mod\)。
   4. 計算目標函數（平方誤差和）。
   5. 顯示擬合結果和誤差。

5. Riedel 方程：非線性最小二乘法
   ```matlab
   %% Riedel Equation: Nonlinear Least Squares
   riedel_fun = @(params, T) 10.^(params(1) + params(2)./T + params(3)*log10(T) + params(4)*T.^2);
   p0 = [0, 0, 0, 0]; % Initial guess for A, B, C, and D

   % Perform nonlinear least squares
   [p_riedel, resnorm_riedel, residual_riedel] = lsqcurvefit(riedel_fun, p0, T, log10(P));

   % Extract parameters
   A4 = p_riedel(1);
   B4 = p_riedel(2);
   C4 = p_riedel(3);
   D4 = p_riedel(4);
   J4 = resnorm_riedel;

   % Display results
   fprintf('Nonlinear Least Squares for Riedel Equation:\n');
   fprintf('A = %.4f\n', A4);
   fprintf('B = %.4f\n', B4);
   fprintf('C = %.4f\n', C4);
   fprintf('D = %.4f\n', D4);
   fprintf('Sum of Squared Errors (J4) = %.4f\n\n', J4);
   ```
   這部分程式碼使用非線性最小二乘法來擬合 Riedel 方程：
   1. 定義非線性模型 `riedel_fun`。
   2. 設置初始參數猜測 \(p0\)。
   3. 使用 `lsqcurvefit` 函數進行擬合。
   4. 提取擬合後的參數 \(A4\)、\(B4\)、\(C4\) 和 \(D4\)。
   5. 計算並顯示誤差。


### Result

```
Linear Least Squares for Clapeyron Equation:
A = 8.7520
B = -2035.3310
Sum of Squared Errors (J1) = 50307.5085
Local minimum found.

Nonlinear Least Squares for Clapeyron Equation:
A = 1.8677
B = -484.2378
Sum of Squared Errors (J2) = 0.7310

Linear Least Squares for Riedel Equation:
A = 216.6854
B = -9317.2030
C = -75.7355
D = 0.0000
Sum of Squared Errors (J3) = 141.4212

Nonlinear Least Squares for Riedel Equation:
A = 595.4824
B = -20462.2702
C = -218.3222
D = 0.0002
Sum of Squared Errors (J4) = 0.0649
```