# 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擬合並不精確:

## 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

## 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
```