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