% --- Comparison of Root-Finding Methods V1.2 Functional--- % clc; clear; close all; % % % Define ranges % temp_range = 100:200:3000; % Temperature range (K) % press_range = 20:10:120; % Pressure range (bar) % feq_range = 0.8:0.05:1.3; % Fuel equivalency ratio range % % % Gibbs free energy table % gibbseng=[0 -393.151 -238.921 -113.805 % 100 -393.683 -236.584 -120.239 % 200 -394.085 -232.766 -128.526 % 298.15 -394.389 -228.582 -137.163 % 300 -394.394 -228.5 -137.328 % 400 -394.675 -223.901 -146.338 % 500 -394.939 -219.051 -155.414 % 600 -395.182 -214.007 -164.486 % 700 -395.398 -208.812 -173.518 % 800 -395.586 -203.496 -182.497 % 900 -395.748 -198.083 -191.416 % 1000 -395.886 -192.59 -200.275 % 1100 -396.001 -187.033 -209.075 % 1200 -396.098 -181.425 -217.819 % 1300 -396.177 -175.774 -226.509 % 1400 -396.24 -170.089 -235.149 % 1500 -396.288 -164.376 -243.74 % 1600 -396.323 -158.639 -252.284 % 1700 -396.344 -152.883 -260.784 % 1800 -396.353 -147.111 -269.242 % 1900 -396.349 -141.325 -277.658 % 2000 -396.333 -135.528 -286.034 % 2100 -396.304 -129.721 -294.372 % 2200 -396.262 -123.905 -302.672 % 2300 -396.209 -118.082 -310.936 % 2400 -396.142 -112.252 -319.165 % 2500 -396.062 -106.416 -327.358 % 2600 -395.969 -100.575 -335.517 % 2700 -395.862 -94.729 -343.643 % 2800 -395.742 -88.878 -351.736 % 2900 -395.609 -83.023 -359.797 % 3000 -395.461 -77.163 -367.826 % 3100 -395.298 -71.298 -375.824 % 3200 -395.122 -65.43 -383.792 % 3300 -394.932 -59.558 -391.73 % 3400 -394.726 -53.681 -399.638 % 3500 -394.506 -47.801 -407.517 % 3600 -394.271 -41.916 -415.366 % 3700 -394.022 -36.027 -423.188 % 3800 -393.756 -30.133 -430.981 % 3900 -393.477 -24.236 -438.746 % 4000 -393.183 -18.334 -446.485 % 4100 -392.874 -12.427 -454.196 % 4200 -392.55 -6.516 -461.881 % 4300 -392.21 -0.6 -469.539 % 4400 -391.857 5.32 -477.171 % 4500 -391.488 11.245 -484.778 % 4600 -391.105 17.175 -492.359 % 4700 -390.706 23.111 -499.915 % 4800 -390.292 29.052 -507.446 % 4900 -389.862 34.998 -514.951 % 5000 -389.419 40.949 -522.433 % 5100 -388.959 46.906 -529.89 % 5200 -388.486 52.869 -537.323 % 5300 -387.996 58.838 -544.732 % 5400 -387.493 64.811 -552.118 % 5500 -386.974 70.791 -559.48 % 5600 -386.439 76.777 -566.818 % 5700 -385.89 82.769 -574.133 % 5800 -385.324 88.767 -581.424 % 5900 -384.745 94.77 -588.694 % 6000 -384.148 100.78 -595.939]; % % % Randomly select values for comparison V1 Functional % % random_temp = temp_range(randi(length(temp_range))); % Random temperature % % random_press = press_range(randi(length(press_range))); % Random pressure % % random_feq = feq_range(randi(length(feq_range))); % Random equivalency ratio % % % % fprintf('\nRandomly selected values for comparison:\n'); % % fprintf(' Temperature: %.2f K\n', random_temp); % % fprintf(' Pressure: %.2f bar\n', random_press); % % fprintf(' Fuel Equiv. Ratio: %.2f\n', random_feq); % % % % % --------------------------- % % % Bisection Method % % % --------------------------- % % fprintf('\n--- Using Bisection Method ---\n'); % % xLow = 0.1; % % xHigh = 1; % % tol = 1e-12; % % maxIter = 100; % % iterCount_bisection = 0; % % % % % Start timing % % tic; % % fLow = F_of_x(xLow, gibbseng, random_temp, random_press, random_feq); % % fHigh = F_of_x(xHigh, gibbseng, random_temp, random_press, random_feq); % % % % % Expand bounds if necessary % % while fLow * fHigh > 0 % % xHigh = xHigh + 0.5; % % fHigh = F_of_x(xHigh, gibbseng, random_temp, random_press, random_feq); % % if xHigh > 5 % % error('Bisection failed to bracket the root.'); % % end % % end % % % % % Perform Bisection % % xMid_bisection = 0; % % fMid_bisection = 0; % % while iterCount_bisection < maxIter % % iterCount_bisection = iterCount_bisection + 1; % % xMid_bisection = (xLow + xHigh) / 2; % % fMid_bisection = F_of_x(xMid_bisection, gibbseng, random_temp, random_press, random_feq); % % % % % Check for convergence % % if abs(fMid_bisection) < tol % % break; % % end % % % % % Update bounds % % if fLow * fMid_bisection < 0 % % xHigh = xMid_bisection; % % fHigh = fMid_bisection; % % else % % xLow = xMid_bisection; % % fLow = fMid_bisection; % % end % % end % % bisection_time = toc; % % % % % Calculate Oxygen Balance for Bisection % % oxygen_balance_bisection = oxygenBalance(xMid_bisection, gibbseng, random_temp, random_press, random_feq); % % % % % % % --------------------------- % % % fzero Method % % % --------------------------- % % fprintf('\n--- Using fzero ---\n'); % % options = optimset('TolFun', 1e-12, 'Display', 'off'); % % tic; % % [x_fzero, fval_fzero, exitflag_fzero, output_fzero] = ... % % fzero(@(x) F_of_x(x, gibbseng, random_temp, random_press, random_feq), [0.1, 1], options); % % fzero_time = toc; % % % % % --------------------------- % % % Original Method % % % --------------------------- % % fprintf('\n--- Using Original Function ---\n'); % % tic; % % [n5ntVal_orig, errorVal_orig, iterCount_orig] = mydissfuniii_original( ... % % random_temp, random_press, random_feq, gibbseng); % Updated original function % % original_time = toc; % % % % % --------------------------- % % % Create Comparison Table % % % --------------------------- % % method_data = { % % 'Bisection', random_temp, random_press, random_feq, xMid_bisection, abs(fMid_bisection), iterCount_bisection, bisection_time; % % 'fzero', random_temp, random_press, random_feq, x_fzero, abs(fval_fzero), output_fzero.iterations, fzero_time; % % 'Original Function', random_temp, random_press, random_feq, n5ntVal_orig, abs(errorVal_orig), iterCount_orig, original_time % % }; % % % % comparison_table = cell2table(method_data, ... % % 'VariableNames', {'Method', 'Temperature', 'Pressure', 'Fuel Eq. Ratio', 'Root', 'Error', 'Total Iterations', 'Comp Time (s)'}); % % % % disp(comparison_table); % % % % % --------------------------- % % % Nested Functions % % % --------------------------- % % % % function FF = F_of_x(x, gibbseng, temperature, pressure, fuelequiv) % % % Interpolate Gibbs free energy (J/mol) % % gco2 = interp1(gibbseng(:,1), gibbseng(:,2), temperature) * 1e3; % % gh2o = interp1(gibbseng(:,1), gibbseng(:,3), temperature) * 1e3; % % gco = interp1(gibbseng(:,1), gibbseng(:,4), temperature) * 1e3; % % % % % Equilibrium constants % % kco2 = exp(-(gco2 - gco) * 1e3 / (8314.46 * temperature)); % % kh2o = exp(-(gh2o) * 1e3 / (8314.46 * temperature)); % % % % % Compute F(x) % % a = sqrt(x * pressure); % % n3 = fuelequiv / (kco2 * a + 1); % % n4 = fuelequiv / (kh2o * a + 1); % % n5 = (1 - 2*n3 - n4) / 2; % % nt = n3 + n4 + n5 + 3.76; % Total moles (simplified) % % FF = x - (n5 / nt); % % end % % % % function [n5ntVal, errorVal, iterCount] = mydissfuniii_original(temp, press, feq, gibbseng) % % % Initialize % % n5ntvalue = 0.1; % % errorVal = 1; % % iterCount = 0; % % % % % Interpolate Gibbs free energy (J/mol) % % gco2 = interp1(gibbseng(:,1), gibbseng(:,2), temp) * 1e3; % % gh2o = interp1(gibbseng(:,1), gibbseng(:,3), temp) * 1e3; % % gco = interp1(gibbseng(:,1), gibbseng(:,4), temp) * 1e3; % % % % % Equilibrium constants % % kco2 = exp(-(gco2 - gco) * 1e3 / (8314.46 * temp)); % % kh2o = exp(-(gh2o) * 1e3 / (8314.46 * temp)); % % % % % Iterative process % % while errorVal > 0 % % iterCount = iterCount + 1; % % a = sqrt(n5ntvalue * press); % % n3 = feq / (kco2 * a + 1); % % n4 = feq / (kh2o * a + 1); % % n5 = (1 - 2*n3 - n4) / 2; % % nt = n3 + n4 + n5 + 3.76; % % errorVal = n5ntvalue - (n5 / nt); % % n5ntvalue = n5ntvalue - 0.001; % % if n5ntvalue < 0 % % break; % % end % % end % % % % % Output results % % n5ntVal = n5ntvalue; % % end % % % % % % % Randomly select values for comparison V2 Functional % % random_temp = temp_range(randi(length(temp_range))); % Random temperature % % random_press = press_range(randi(length(press_range))); % Random pressure % % random_feq = feq_range(randi(length(feq_range))); % Random equivalency ratio % % % % fprintf('\nRandomly selected values for comparison:\n'); % % fprintf(' Temperature: %.2f K\n', random_temp); % % fprintf(' Pressure: %.2f bar\n', random_press); % % fprintf(' Fuel Equiv. Ratio: %.2f\n', random_feq); % % % % % --------------------------- % % % Bisection Method % % % --------------------------- % % fprintf('\n--- Using Bisection Method ---\n'); % % xLow = 0.1; % % xHigh = 4; % % tol = 1e-9; % % maxIter = 100; % % iterCount_bisection = 0; % % % % % Start timing % % tic; % % fLow = F_of_x(xLow, gibbseng, random_temp, random_press, random_feq); % % fHigh = F_of_x(xHigh, gibbseng, random_temp, random_press, random_feq); % % % % % Expand bounds if necessary % % while fLow * fHigh > 0 % % xHigh = xHigh + 0.5; % % fHigh = F_of_x(xHigh, gibbseng, random_temp, random_press, random_feq); % % if xHigh > 5 % % error('Bisection failed to bracket the root.'); % % end % % end % % % % % Perform Bisection % % xMid_bisection = 0; % % fMid_bisection = 0; % % while iterCount_bisection < maxIter % % iterCount_bisection = iterCount_bisection + 1; % % xMid_bisection = (xLow + xHigh) / 2; % % fMid_bisection = F_of_x(xMid_bisection, gibbseng, random_temp, random_press, random_feq); % % % % % Check for convergence % % if abs(fMid_bisection) < tol % % break; % % end % % % % % Update bounds % % if fLow * fMid_bisection < 0 % % xHigh = xMid_bisection; % % fHigh = fMid_bisection; % % else % % xLow = xMid_bisection; % % fLow = fMid_bisection; % % end % % end % % bisection_time = toc; % % % % % Calculate Oxygen Balance for Bisection % % oxygen_balance_bisection = oxygenBalance(xMid_bisection, gibbseng, random_temp, random_press, random_feq); % % % % % --------------------------- % % % fzero Method % % % --------------------------- % % fprintf('\n--- Using fzero ---\n'); % % options = optimset('TolFun', 1e-9, 'Display', 'off'); % % tic; % % [x_fzero, fval_fzero, exitflag_fzero, output_fzero] = ... % % fzero(@(x) F_of_x(x, gibbseng, random_temp, random_press, random_feq), [0.1, 4], options); % % fzero_time = toc; % % % % % Calculate Oxygen Balance for fzero % % oxygen_balance_fzero = oxygenBalance(x_fzero, gibbseng, random_temp, random_press, random_feq); % % % % % --------------------------- % % % Original Method % % % --------------------------- % % fprintf('\n--- Using Original Function ---\n'); % % tic; % % [n5ntVal_orig, errorVal_orig, iterCount_orig] = mydissfuniii_original( ... % % random_temp, random_press, random_feq, gibbseng); % Updated original function % % original_time = toc; % % % % % Calculate Oxygen Balance for Original % % oxygen_balance_original = oxygenBalance(n5ntVal_orig, gibbseng, random_temp, random_press, random_feq); % % % % % --------------------------- % % % Create Comparison Table % % % --------------------------- % % method_data = { % % 'Bisection', random_temp, random_press, random_feq, xMid_bisection, abs(fMid_bisection), iterCount_bisection, bisection_time, oxygen_balance_bisection; % % 'fzero', random_temp, random_press, random_feq, x_fzero, abs(fval_fzero), output_fzero.iterations, fzero_time, oxygen_balance_fzero; % % 'Original Function', random_temp, random_press, random_feq, n5ntVal_orig, abs(errorVal_orig), iterCount_orig, original_time, oxygen_balance_original % % }; % % % % comparison_table = cell2table(method_data, ... % % 'VariableNames', {'Method', 'Temperature', 'Pressure', 'Fuel Eq. Ratio', 'Root', 'Error', 'Total Iterations', 'Comp Time (s)', 'Oxygen Balance'}); % % % % % Display the table % % disp(comparison_table); % % % % % --------------------------- % % % Nested Functions % % % --------------------------- % % % % function FF = F_of_x(x, gibbseng, temperature, pressure, fuelequiv) % % % Interpolate Gibbs free energy (J/mol) % % gco2 = interp1(gibbseng(:,1), gibbseng(:,2), temperature) * 1e3; % % gh2o = interp1(gibbseng(:,1), gibbseng(:,3), temperature) * 1e3; % % gco = interp1(gibbseng(:,1), gibbseng(:,4), temperature) * 1e3; % % % % % Equilibrium constants % % kco2 = exp(-(gco2 - gco) * 1e3 / (8314.46 * temperature)); % % kh2o = exp(-(gh2o) * 1e3 / (8314.46 * temperature)); % % % % % Compute F(x) % % a = sqrt(x * pressure); % % n3 = fuelequiv / (kco2 * a + 1); % % n4 = fuelequiv / (kh2o * a + 1); % % n5 = (1 - 2*n3 - n4) / 2; % % nt = n3 + n4 + n5 + 3.76; % Total moles (simplified) % % FF = x - (n5 / nt); % % end % % % % function [n5ntVal, errorVal, iterCount] = mydissfuniii_original(temp, press, feq, gibbseng) % % % Initialize % % n5ntvalue = 0.1; % % errorVal = 1; % % iterCount = 0; % % % % % Interpolate Gibbs free energy (J/mol) % % gco2 = interp1(gibbseng(:,1), gibbseng(:,2), temp) * 1e3; % % gh2o = interp1(gibbseng(:,1), gibbseng(:,3), temp) * 1e3; % % gco = interp1(gibbseng(:,1), gibbseng(:,4), temp) * 1e3; % % % % % Equilibrium constants % % kco2 = exp(-(gco2 - gco) * 1e3 / (8314.46 * temp)); % % kh2o = exp(-(gh2o) * 1e3 / (8314.46 * temp)); % % % % % Iterative process % % while errorVal > 0 % % iterCount = iterCount + 1; % % a = sqrt(n5ntvalue * press); % % n3 = feq / (kco2 * a + 1); % % n4 = feq / (kh2o * a + 1); % % n5 = (1 - 2*n3 - n4) / 2; % % nt = n3 + n4 + n5 + 3.76; % % errorVal = n5ntvalue - (n5 / nt); % % n5ntvalue = n5ntvalue - 0.001; % % if n5ntvalue < 0 % % break; % % end % % end % % % % % Output results % % n5ntVal = n5ntvalue; % % end % % % % function oxygen_balance = oxygenBalance(root, gibbseng, temp, press, feq) % % % Compute oxygen balance for given root % % a = sqrt(root * press); % % gco2 = interp1(gibbseng(:,1), gibbseng(:,2), temp) * 1e3; % % gco = interp1(gibbseng(:,1), gibbseng(:,4), temp) * 1e3; % % kco2 = exp(-(gco2 - gco) * 1e3 / (8314.46 * temp)); % % oxygen_balance = (feq / (kco2 * a + 1)) - root; % Simplified example % % end % % % Randomly select values for comparison % random_temp = temp_range(randi(length(temp_range))); % Random temperature % random_press = press_range(randi(length(press_range))); % Random pressure % random_feq = feq_range(randi(length(feq_range))); % Random equivalency ratio % % fprintf('\nRandomly selected values for comparison:\n'); % fprintf(' Temperature: %.2f K\n', random_temp); % fprintf(' Pressure: %.2f bar\n', random_press); % fprintf(' Fuel Equiv. Ratio: %.2f\n', random_feq); % % % --------------------------- % % Bisection Method % % --------------------------- % fprintf('\n--- Using Bisection Method ---\n'); % xLow = 0; % xHigh = 0.5; % tol = 1e-9; % maxIter = 100; % iterCount_bisection = 0; % % % Start timing % tic; % fLow = F_of_x(xLow, gibbseng, random_temp, random_press, random_feq); % fHigh = F_of_x(xHigh, gibbseng, random_temp, random_press, random_feq); % % % Expand bounds if necessary % while fLow * fHigh > 0 % xHigh = xHigh + 0.5; % fHigh = F_of_x(xHigh, gibbseng, random_temp, random_press, random_feq); % if xHigh > 5 % error('Bisection failed to bracket the root.'); % end % end % % % Perform Bisection % xMid_bisection = 0; % fMid_bisection = 0; % while iterCount_bisection < maxIter % iterCount_bisection = iterCount_bisection + 1; % xMid_bisection = (xLow + xHigh) / 2; % fMid_bisection = F_of_x(xMid_bisection, gibbseng, random_temp, random_press, random_feq); % % % Check for convergence % if abs(fMid_bisection) < tol % break; % end % % % Update bounds % if fLow * fMid_bisection < 0 % xHigh = xMid_bisection; % fHigh = fMid_bisection; % else % xLow = xMid_bisection; % fLow = fMid_bisection; % end % end % bisection_time = toc; % % % Calculate Oxygen Balance for Bisection % oxygen_balance_bisection = calculate_oxygen_balance(xMid_bisection, random_temp, random_press, random_feq, gibbseng); % % % --------------------------- % % fzero Method % % --------------------------- % fprintf('\n--- Using fzero ---\n'); % options = optimset('TolFun', 1e-9, 'Display', 'off'); % tic; % [x_fzero, fval_fzero, exitflag_fzero, output_fzero] = ... % fzero(@(x) F_of_x(x, gibbseng, random_temp, random_press, random_feq), [0, 0.5], options); % fzero_time = toc; % % % Calculate Oxygen Balance for fzero % oxygen_balance_fzero = calculate_oxygen_balance(x_fzero, random_temp, random_press, random_feq, gibbseng); % % % --------------------------- % % Original Method % % --------------------------- % fprintf('\n--- Using Original Function ---\n'); % tic; % [n5ntVal_orig, errorVal_orig, iterCount_orig] = mydissfuniii_original( ... % random_temp, random_press, random_feq, gibbseng); % Updated original function % original_time = toc; % % % Calculate Oxygen Balance for Original % oxygen_balance_original = calculate_oxygen_balance(n5ntVal_orig, random_temp, random_press, random_feq, gibbseng); % % % --------------------------- % % Create Comparison Table % % --------------------------- % method_data = { % 'Bisection', random_temp, random_press, random_feq, xMid_bisection, abs(fMid_bisection), iterCount_bisection, bisection_time, oxygen_balance_bisection; % 'fzero', random_temp, random_press, random_feq, x_fzero, abs(fval_fzero), output_fzero.iterations, fzero_time, oxygen_balance_fzero; % 'Original Function', random_temp, random_press, random_feq, n5ntVal_orig, abs(errorVal_orig), iterCount_orig, original_time, oxygen_balance_original % }; % % comparison_table = cell2table(method_data, ... % 'VariableNames', {'Method', 'Temperature', 'Pressure', 'Fuel Eq. Ratio', 'Root', 'Error', 'Total Iterations', 'Comp Time (s)', 'Oxygen Balance'}); % % % Display the table % disp(comparison_table); % % % --------------------------- % % Nested Functions % % --------------------------- % % function FF = F_of_x(x, gibbseng, temperature, pressure, fuelequiv) % % Interpolate Gibbs free energy (J/mol) % gco2 = interp1(gibbseng(:,1), gibbseng(:,2), temperature) * 1e3; % gh2o = interp1(gibbseng(:,1), gibbseng(:,3), temperature) * 1e3; % gco = interp1(gibbseng(:,1), gibbseng(:,4), temperature) * 1e3; % % % Equilibrium constants % kco2 = exp(-(gco2 - gco) * 1e3 / (8314.46 * temperature)); % kh2o = exp(-(gh2o) * 1e3 / (8314.46 * temperature)); % % % Compute F(x) % a = sqrt(x * pressure); % n3 = fuelequiv / (kco2 * a + 1); % n4 = fuelequiv / (kh2o * a + 1); % n5 = max(0, (1 - 2*n3 - n4) / 2); % Ensure n5 is non-negative % nt = n3 + n4 + n5 + 3.76; % Total moles (simplified) % FF = x - (n5 / nt); % end % % function oxygen_balance = calculate_oxygen_balance(root, temperature, pressure, fuelequiv, gibbseng) % % Compute oxygen balance for given root % a = sqrt(root * pressure); % % % Interpolate Gibbs free energy (J/mol) % gco2 = interp1(gibbseng(:,1), gibbseng(:,2), temperature) * 1e3; % gh2o = interp1(gibbseng(:,1), gibbseng(:,3), temperature) * 1e3; % gco = interp1(gibbseng(:,1), gibbseng(:,4), temperature) * 1e3; % % % Equilibrium constants % kco2 = exp(-(gco2 - gco) * 1e3 / (8314.46 * temperature)); % % % Oxygen balance equation % n3 = fuelequiv / (kco2 * a + 1); % oxygen_balance = max(0, 1 - 2 * n3); % Ensure non-negative balance % end % % --- Comparison of Root-Finding Methods V2 Functional--- clc; clear; close all; % Define ranges temp_range = 100:200:3000; % Temperature range (K) press_range = 20:10:120; % Pressure range (bar) feq_range = 0.8:0.05:1.3; % Fuel equivalency ratio range % Randomly select values for comparison random_temp = temp_range(randi(length(temp_range))); % Random temperature random_press = press_range(randi(length(press_range))); % Random pressure random_feq = feq_range(randi(length(feq_range))); % Random equivalency ratio fprintf('\nRandomly selected values for comparison:\n'); fprintf(' Temperature: %.2f K\n', random_temp); fprintf(' Pressure: %.2f bar\n', random_press); fprintf(' Fuel Equiv. Ratio: %.2f\n', random_feq); % --------------------------- % Bisection Method % --------------------------- fprintf('\n--- Using Bisection Method ---\n'); xLow = 0.1; xHigh = 4; tol = 1e-9; maxIter = 100; iterCount_bisection = 0; % Start timing tic; fLow = F_of_x(xLow, random_temp, random_press, random_feq); fHigh = F_of_x(xHigh, random_temp, random_press, random_feq); % Expand bounds if necessary while fLow * fHigh > 0 xHigh = xHigh + 0.5; fHigh = F_of_x(xHigh, random_temp, random_press, random_feq); if xHigh > 5 error('Bisection failed to bracket the root.'); end end % Perform Bisection xMid_bisection = 0; fMid_bisection = 0; while iterCount_bisection < maxIter iterCount_bisection = iterCount_bisection + 1; xMid_bisection = (xLow + xHigh) / 2; fMid_bisection = F_of_x(xMid_bisection, random_temp, random_press, random_feq); % Check for convergence if abs(fMid_bisection) < tol break; end % Update bounds if fLow * fMid_bisection < 0 xHigh = xMid_bisection; fHigh = fMid_bisection; else xLow = xMid_bisection; fLow = fMid_bisection; end end bisection_time = toc; % Calculate Oxygen Balance for Bisection oxygen_balance_bisection = calculate_oxygen_balance(xMid_bisection, random_temp, random_press, random_feq); % --------------------------- % fzero Method % --------------------------- fprintf('\n--- Using fzero ---\n'); options = optimset('TolFun', 1e-9, 'Display', 'off'); tic; [x_fzero, fval_fzero, exitflag_fzero, output_fzero] = ... fzero(@(x) F_of_x(x, random_temp, random_press, random_feq), [0.1, 4], options); fzero_time = toc; % Calculate Oxygen Balance for fzero oxygen_balance_fzero = calculate_oxygen_balance(x_fzero, random_temp, random_press, random_feq); % --------------------------- % Original Method (mydissfuniii) % --------------------------- fprintf('\n--- Using Original Function ---\n'); tic; [out_orig, error_orig, iterCount_orig] = mydissfuniii(random_temp, random_press, random_feq); % Modified mydissfuniii original_time = toc; % Calculate Oxygen Balance for Original oxygen_balance_original = calculate_oxygen_balance(out_orig, random_temp, random_press, random_feq); % --------------------------- % Create Comparison Table % --------------------------- method_data = { 'Bisection', random_temp, random_press, random_feq, xMid_bisection, abs(fMid_bisection), iterCount_bisection, bisection_time, oxygen_balance_bisection; 'fzero', random_temp, random_press, random_feq, x_fzero, abs(fval_fzero), output_fzero.iterations, fzero_time, oxygen_balance_fzero; 'Original Function', random_temp, random_press, random_feq, out_orig, abs(error_orig), iterCount_orig, original_time, oxygen_balance_original }; comparison_table = cell2table(method_data, ... 'VariableNames', {'Method', 'Temperature', 'Pressure', 'Fuel Eq. Ratio', 'Root', 'Error', 'Total Iterations', 'Comp Time (s)', 'Oxygen Balance'}); % Display the table disp(comparison_table); % --------------------------- % Nested Functions % --------------------------- function FF = F_of_x(x, temperature, pressure, fuelequiv) % Interpolate Gibbs free energy gibbseng = load_gibbs(); % Load Gibbs data gco2 = interp1(gibbseng(:,1), gibbseng(:,2), temperature) * 1e3; gh2o = interp1(gibbseng(:,1), gibbseng(:,3), temperature) * 1e3; gco = interp1(gibbseng(:,1), gibbseng(:,4), temperature) * 1e3; % Equilibrium constants kco2 = exp(-(gco2 - gco) * 1e3 / (8314.46 * temperature)); kh2o = exp(-(gh2o) * 1e3 / (8314.46 * temperature)); % Compute F(x) a = sqrt(x * pressure); n3 = fuelequiv / (kco2 * a + 1); n4 = fuelequiv / (kh2o * a + 1); n5 = max(0, (1 - 2*n3 - n4) / 2); % Ensure n5 is non-negative nt = n3 + n4 + n5 + 3.76; % Total moles (simplified) FF = x - (n5 / nt); end function [out, errorVal, iterCount] = mydissfuniii(temperature, pressure, fuelequiv) % Original mydissfuniii function with added iteration tracking iterCount = 0; n5ntvalue = 0.1; errorVal = 1; % Interpolate Gibbs free energy gibbseng = load_gibbs(); gco2 = interp1(gibbseng(:,1), gibbseng(:,2), temperature) * 1e3; gh2o = interp1(gibbseng(:,1), gibbseng(:,3), temperature) * 1e3; gco = interp1(gibbseng(:,1), gibbseng(:,4), temperature) * 1e3; % Equilibrium constants kco2 = exp(-(gco2 - gco) * 1e3 / (8314.46 * temperature)); kh2o = exp(-(gh2o) * 1e3 / (8314.46 * temperature)); % Iterative process while errorVal > 0 iterCount = iterCount + 1; a = sqrt(n5ntvalue * pressure); n3 = fuelequiv / (kco2 * a + 1); n4 = fuelequiv / (kh2o * a + 1); n5 = (1 - 2*n3 - n4) / 2; nt = n3 + n4 + n5 + 3.76; errorVal = n5ntvalue - (n5 / nt); n5ntvalue = n5ntvalue - 0.001; if n5ntvalue < 0 break; end end out = n5ntvalue; end function oxygen_balance = calculate_oxygen_balance(root, temperature, pressure, fuelequiv) % Compute oxygen balance for given root a = sqrt(root * pressure); gibbseng = load_gibbs(); gco2 = interp1(gibbseng(:,1), gibbseng(:,2), temperature) * 1e3; gco = interp1(gibbseng(:,1), gibbseng(:,4), temperature) * 1e3; kco2 = exp(-(gco2 - gco) * 1e3 / (8314.46 * temperature)); % Oxygen balance equation n3 = fuelequiv / (kco2 * a + 1); oxygen_balance = max(0, 1 - 2 * n3); % Ensure non-negative balance end function gibbseng = load_gibbs() gibbseng=[0 -393.151 -238.921 -113.805 100 -393.683 -236.584 -120.239 200 -394.085 -232.766 -128.526 298.15 -394.389 -228.582 -137.163 300 -394.394 -228.5 -137.328 400 -394.675 -223.901 -146.338 500 -394.939 -219.051 -155.414 600 -395.182 -214.007 -164.486 700 -395.398 -208.812 -173.518 800 -395.586 -203.496 -182.497 900 -395.748 -198.083 -191.416 1000 -395.886 -192.59 -200.275 1100 -396.001 -187.033 -209.075 1200 -396.098 -181.425 -217.819 1300 -396.177 -175.774 -226.509 1400 -396.24 -170.089 -235.149 1500 -396.288 -164.376 -243.74 1600 -396.323 -158.639 -252.284 1700 -396.344 -152.883 -260.784 1800 -396.353 -147.111 -269.242 1900 -396.349 -141.325 -277.658 2000 -396.333 -135.528 -286.034 2100 -396.304 -129.721 -294.372 2200 -396.262 -123.905 -302.672 2300 -396.209 -118.082 -310.936 2400 -396.142 -112.252 -319.165 2500 -396.062 -106.416 -327.358 2600 -395.969 -100.575 -335.517 2700 -395.862 -94.729 -343.643 2800 -395.742 -88.878 -351.736 2900 -395.609 -83.023 -359.797 3000 -395.461 -77.163 -367.826 3100 -395.298 -71.298 -375.824 3200 -395.122 -65.43 -383.792 3300 -394.932 -59.558 -391.73 3400 -394.726 -53.681 -399.638 3500 -394.506 -47.801 -407.517 3600 -394.271 -41.916 -415.366 3700 -394.022 -36.027 -423.188 3800 -393.756 -30.133 -430.981 3900 -393.477 -24.236 -438.746 4000 -393.183 -18.334 -446.485 4100 -392.874 -12.427 -454.196 4200 -392.55 -6.516 -461.881 4300 -392.21 -0.6 -469.539 4400 -391.857 5.32 -477.171 4500 -391.488 11.245 -484.778 4600 -391.105 17.175 -492.359 4700 -390.706 23.111 -499.915 4800 -390.292 29.052 -507.446 4900 -389.862 34.998 -514.951 5000 -389.419 40.949 -522.433 5100 -388.959 46.906 -529.89 5200 -388.486 52.869 -537.323 5300 -387.996 58.838 -544.732 5400 -387.493 64.811 -552.118 5500 -386.974 70.791 -559.48 5600 -386.439 76.777 -566.818 5700 -385.89 82.769 -574.133 5800 -385.324 88.767 -581.424 5900 -384.745 94.77 -588.694 6000 -384.148 100.78 -595.939]; end % --- Comparison of Root-Finding Methods for All Cases V3 Not functional--- % clc; clear; close all; % % % Define ranges % temp_range = 100:200:3000; % Temperature range (K) % press_range = 20:10:120; % Pressure range (bar) % feq_range = 0.8:0.05:1.3; % Fuel equivalency ratio range % % % Preallocate results for performance % nTemps = length(temp_range); % nPress = length(press_range); % nFeq = length(feq_range); % nCases = nTemps * nPress * nFeq; % % % Initialize accumulators % iterCounts = zeros(nCases, 3); % Iterations for each method % errors = zeros(nCases, 3); % Errors for each method % times = zeros(1, 3); % Computation times for each method % % % Indices for iterating % caseIndex = 1; % % % --------------------------- % % Bisection Method for All Cases % % --------------------------- % fprintf('\n--- Running Bisection Method for All Cases ---\n'); % tic; % for temp = temp_range % for press = press_range % for feq = feq_range % xLow = 0.01; % xHigh = 1; % tol = 1e-9; % maxIter = 100; % iterCount = 0; % % % Initialize function bounds % fLow = F_of_x(xLow, temp, press, feq); % fHigh = F_of_x(xHigh, temp, press, feq); % % % Expand bounds if necessary % while fLow * fHigh > 0 % xHigh = xHigh + 0.5; % fHigh = F_of_x(xHigh, temp, press, feq); % if xHigh > 5 % break; % Abort if bounds can't bracket % end % end % % % Perform Bisection % xMid = 0; fMid = 0; % while iterCount < maxIter % iterCount = iterCount + 1; % xMid = (xLow + xHigh) / 2; % fMid = F_of_x(xMid, temp, press, feq); % % % Check for convergence % if abs(fMid) < tol % break; % end % % % Update bounds % if fLow * fMid < 0 % xHigh = xMid; % else % xLow = xMid; % end % end % % % Store results % iterCounts(caseIndex, 1) = iterCount; % errors(caseIndex, 1) = abs(fMid); % caseIndex = caseIndex + 1; % end % end % end % times(1) = toc; % % % --------------------------- % % fzero Method % % --------------------------- % fprintf('\n--- Using fzero ---\n'); % % % Define the initial interval % xLow = 0.01; % xHigh = 1; % options = optimset('TolFun', 1e-9, 'Display', 'off'); % % % Check the signs at the endpoints of the interval % fLow = F_of_x(xLow, temp, press, feq); % fHigh = F_of_x(xHigh, temp, press, feq); % % % Expand the interval dynamically if the signs are the same % while fLow * fHigh > 0 % xHigh = xHigh + 0.5; % fHigh = F_of_x(xHigh, temp, press, feq); % % % Stop expanding if no root is found % if xHigh > 5 % fprintf('Unable to find a valid interval for fzero.\n'); % root_fzero = NaN; % fval_fzero = NaN; % iterCount_fzero = NaN; % fzero_time = NaN; % return; % end % end % % % If a valid interval is found, proceed with fzero % if fLow * fHigh < 0 % tic; % [root_fzero, fval_fzero, exitflag_fzero, output_fzero] = ... % fzero(@(x) F_of_x(x, temp, press, feq), [xLow, xHigh], options); % fzero_time = toc; % % % Retrieve the number of iterations % if isfield(output_fzero, 'iterations') && ~isempty(output_fzero.iterations) % iterCount_fzero = output_fzero.iterations; % else % iterCount_fzero = NaN; % Handle cases where iterations are not recorded % end % else % % If no valid interval is found, mark as failed % root_fzero = NaN; % fval_fzero = NaN; % iterCount_fzero = NaN; % fzero_time = NaN; % end % % % Debugging Output % fprintf('fzero Results:\n'); % fprintf(' Root: %.6f\n', root_fzero); % fprintf(' Error: %.2e\n', fval_fzero); % fprintf(' Iterations: %d\n', iterCount_fzero); % fprintf(' Computation Time: %.6f s\n', fzero_time); % % % % --------------------------- % % Original Method for All Cases % % --------------------------- % fprintf('\n--- Running Original Function for All Cases ---\n'); % caseIndex = 1; % Reset index % tic; % for temp = temp_range % for press = press_range % for feq = feq_range % % Solve with original method % [root, error, iterCount] = mydissfuniii(temp, press, feq); % % % Store results % iterCounts(caseIndex, 3) = iterCount; % errors(caseIndex, 3) = abs(error); % caseIndex = caseIndex + 1; % end % end % end % times(3) = toc; % % % --------------------------- % % Generate Comparison Table % % --------------------------- % method_names = {'Bisection', 'fzero', 'Original Function'}; % total_iters = sum(iterCounts, 1); % Total iterations for each method % avg_errors = mean(errors, 1); % Average error for each method % % comparison_table = table(method_names', times', total_iters', avg_errors', ... % 'VariableNames', {'Method', 'Total_Time_s', 'Total_Iterations', 'Avg_Error'}); % % % Display the table % disp(comparison_table); % % % --------------------------- % % Nested Functions % % --------------------------- % % function FF = F_of_x(x, temperature, pressure, fuelequiv) % % Interpolate Gibbs free energy % gibbseng = load_gibbs(); % Load Gibbs data % gco2 = interp1(gibbseng(:,1), gibbseng(:,2), temperature) * 1e3; % gh2o = interp1(gibbseng(:,1), gibbseng(:,3), temperature) * 1e3; % gco = interp1(gibbseng(:,1), gibbseng(:,4), temperature) * 1e3; % % % Equilibrium constants % kco2 = exp(-(gco2 - gco) * 1e3 / (8314.46 * temperature)); % kh2o = exp(-(gh2o) * 1e3 / (8314.46 * temperature)); % % % Compute F(x) % a = sqrt(x * pressure); % n3 = fuelequiv / (kco2 * a + 1); % n4 = fuelequiv / (kh2o * a + 1); % n5 = max(0, (1 - 2*n3 - n4) / 2); % Ensure n5 is non-negative % nt = n3 + n4 + n5 + 3.76; % Total moles (simplified) % FF = x - (n5 / nt); % end % % function [out, errorVal, iterCount] = mydissfuniii(temperature, pressure, fuelequiv) % % Original mydissfuniii function with added iteration tracking % iterCount = 0; % n5ntvalue = 0.1; % errorVal = 1; % % % Interpolate Gibbs free energy % gibbseng = load_gibbs(); % gco2 = interp1(gibbseng(:,1), gibbseng(:,2), temperature) * 1e3; % gh2o = interp1(gibbseng(:,1), gibbseng(:,3), temperature) * 1e3; % gco = interp1(gibbseng(:,1), gibbseng(:,4), temperature) * 1e3; % % % Equilibrium constants % kco2 = exp(-(gco2 - gco) * 1e3 / (8314.46 * temperature)); % kh2o = exp(-(gh2o) * 1e3 / (8314.46 * temperature)); % % % Iterative process % while errorVal > 0 % iterCount = iterCount + 1; % a = sqrt(n5ntvalue * pressure); % n3 = fuelequiv / (kco2 * a + 1); % n4 = fuelequiv / (kh2o * a + 1); % n5 = (1 - 2*n3 - n4) / 2; % nt = n3 + n4 + n5 + 3.76; % errorVal = n5ntvalue - (n5 / nt); % n5ntvalue = n5ntvalue - 0.001; % if n5ntvalue < 0 % break; % end % end % % out = n5ntvalue; % end % % function gibbseng = load_gibbs() % gibbseng=[0 -393.151 -238.921 -113.805 % 100 -393.683 -236.584 -120.239 % 200 -394.085 -232.766 -128.526 % 298.15 -394.389 -228.582 -137.163 % 300 -394.394 -228.5 -137.328 % 400 -394.675 -223.901 -146.338 % 500 -394.939 -219.051 -155.414 % 600 -395.182 -214.007 -164.486 % 700 -395.398 -208.812 -173.518 % 800 -395.586 -203.496 -182.497 % 900 -395.748 -198.083 -191.416 % 1000 -395.886 -192.59 -200.275 % 1100 -396.001 -187.033 -209.075 % 1200 -396.098 -181.425 -217.819 % 1300 -396.177 -175.774 -226.509 % 1400 -396.24 -170.089 -235.149 % 1500 -396.288 -164.376 -243.74 % 1600 -396.323 -158.639 -252.284 % 1700 -396.344 -152.883 -260.784 % 1800 -396.353 -147.111 -269.242 % 1900 -396.349 -141.325 -277.658 % 2000 -396.333 -135.528 -286.034 % 2100 -396.304 -129.721 -294.372 % 2200 -396.262 -123.905 -302.672 % 2300 -396.209 -118.082 -310.936 % 2400 -396.142 -112.252 -319.165 % 2500 -396.062 -106.416 -327.358 % 2600 -395.969 -100.575 -335.517 % 2700 -395.862 -94.729 -343.643 % 2800 -395.742 -88.878 -351.736 % 2900 -395.609 -83.023 -359.797 % 3000 -395.461 -77.163 -367.826 % 3100 -395.298 -71.298 -375.824 % 3200 -395.122 -65.43 -383.792 % 3300 -394.932 -59.558 -391.73 % 3400 -394.726 -53.681 -399.638 % 3500 -394.506 -47.801 -407.517 % 3600 -394.271 -41.916 -415.366 % 3700 -394.022 -36.027 -423.188 % 3800 -393.756 -30.133 -430.981 % 3900 -393.477 -24.236 -438.746 % 4000 -393.183 -18.334 -446.485 % 4100 -392.874 -12.427 -454.196 % 4200 -392.55 -6.516 -461.881 % 4300 -392.21 -0.6 -469.539 % 4400 -391.857 5.32 -477.171 % 4500 -391.488 11.245 -484.778 % 4600 -391.105 17.175 -492.359 % 4700 -390.706 23.111 -499.915 % 4800 -390.292 29.052 -507.446 % 4900 -389.862 34.998 -514.951 % 5000 -389.419 40.949 -522.433 % 5100 -388.959 46.906 -529.89 % 5200 -388.486 52.869 -537.323 % 5300 -387.996 58.838 -544.732 % 5400 -387.493 64.811 -552.118 % 5500 -386.974 70.791 -559.48 % 5600 -386.439 76.777 -566.818 % 5700 -385.89 82.769 -574.133 % 5800 -385.324 88.767 -581.424 % 5900 -384.745 94.77 -588.694 % 6000 -384.148 100.78 -595.939]; % end % %% %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% %Root Finding Methods %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% %% solveAllMethods function solveAllMethods %% 1) Define the function and domain clc; clear;close all f = @(x) 10*log10(x) - x; % Our function; domain x>0 df = @(x) (10/(x*log(10))) - 1; % Analytical derivative xLower = 1e-6; % "0" replaced by a small positive xUpper = 4; % Upper bound tol = 1e-7; % ~6 decimal accuracy maxIter = 500; % Maximum iteration limit % We'll create 7 methods, store results in a struct array of length 7. nMethods = 7; methodsData(nMethods) = struct('methodName',[],'root',NaN,'fVal',NaN,... 'iterations',NaN,'timeSec',NaN); %% For plotting f(x) xx = linspace(0,4,300); ff = arrayfun(f, xx); iMethod = 1; % We'll increment this as we fill methodsData %% 2) Method #1: Bisection methodName = 'Bisection'; figure('Name',methodName,'Color','w'); plot(xx, ff,'b-','LineWidth',2); hold on; grid on; xlabel('x'); ylabel('f(x)'); title(methodName); tic; [rootBis, itBis, xHistB] = bisectionMethod(f, xLower, xUpper, tol, maxIter); tBis = toc; fValBis = fCheck(rootBis); % Plot iteration points & final root plotIteration(xx, ff, f, xHistB, rootBis); methodsData(iMethod).methodName = methodName; methodsData(iMethod).root = rootBis; methodsData(iMethod).fVal = fValBis; methodsData(iMethod).iterations = itBis; methodsData(iMethod).timeSec = tBis; iMethod = iMethod + 1; %% 3) Method #2: Secant methodName = 'Secant'; figure('Name',methodName,'Color','w'); plot(xx, ff,'b-','LineWidth',2); hold on; grid on; xlabel('x'); ylabel('f(x)'); title(methodName); tic; [rootSec, itSec, xHistS] = secantMethod(f, xLower, xUpper, tol, maxIter); tSec = toc; fValSec = fCheck(rootSec); plotIteration(xx, ff, f, xHistS, rootSec); methodsData(iMethod).methodName = methodName; methodsData(iMethod).root = rootSec; methodsData(iMethod).fVal = fValSec; methodsData(iMethod).iterations = itSec; methodsData(iMethod).timeSec = tSec; iMethod = iMethod + 1; %% 4) Method #3: Newton–Raphson (double precision) methodName = 'Newton-Raphson'; figure('Name',methodName,'Color','w'); plot(xx, ff,'b-','LineWidth',2); hold on; grid on; xlabel('x'); ylabel('f(x)'); title(methodName); tic; [rootNewt, itNewt, xHistN] = newtonMethod(f, df, xLower, tol, maxIter); tNewt = toc; fValNewt = fCheck(rootNewt); plotIteration(xx, ff, f, xHistN, rootNewt); methodsData(iMethod).methodName = methodName; methodsData(iMethod).root = rootNewt; methodsData(iMethod).fVal = fValNewt; methodsData(iMethod).iterations = itNewt; methodsData(iMethod).timeSec = tNewt; iMethod = iMethod + 1; %% 5) Method #4: False Position methodName = 'False Position'; figure('Name',methodName,'Color','w'); plot(xx, ff,'b-','LineWidth',2); hold on; grid on; xlabel('x'); ylabel('f(x)'); title(methodName); tic; [rootFals, itFals, xHistF] = falsePositionMethod(f, xLower, xUpper, tol, maxIter); tFals = toc; fValFals = fCheck(rootFals); plotIteration(xx, ff, f, xHistF, rootFals); methodsData(iMethod).methodName = methodName; methodsData(iMethod).root = rootFals; methodsData(iMethod).fVal = fValFals; methodsData(iMethod).iterations = itFals; methodsData(iMethod).timeSec = tFals; iMethod = iMethod + 1; %% 6) Method #5: fzero (Built-in) methodName = 'fzero (Built-in)'; figure('Name',methodName,'Color','w'); plot(xx, ff,'b-','LineWidth',2); hold on; grid on; xlabel('x'); ylabel('f(x)'); title(methodName); t0 = tic; try opts = optimset('Display','off','TolX',tol,'TolFun',tol); rootFz = fzero(f, [xLower, xUpper], opts); tFz = toc(t0); fValFz = fCheck(rootFz); plotIteration(xx, ff, f, [xLower, rootFz], rootFz); itFz = 1; % Not tracked by fzero easily methodsData(iMethod).methodName = methodName; methodsData(iMethod).root = rootFz; methodsData(iMethod).fVal = fValFz; methodsData(iMethod).iterations = itFz; methodsData(iMethod).timeSec = tFz; catch tFz = toc(t0); methodsData(iMethod).methodName = methodName; methodsData(iMethod).root = NaN; methodsData(iMethod).fVal = NaN; methodsData(iMethod).iterations = NaN; methodsData(iMethod).timeSec = tFz; end iMethod = iMethod + 1; %% 7) Method #6: fsolve (Built-in) methodName = 'fsolve (Built-in)'; if exist('fsolve','file')==2 figure('Name',methodName,'Color','w'); plot(xx, ff,'b-','LineWidth',2); hold on; grid on; xlabel('x'); ylabel('f(x)'); title(methodName); t0 = tic; try opts2 = optimoptions('fsolve','Display','off','TolX',tol,... 'TolFun',tol,'MaxIterations',maxIter); rootFs = fsolve(f, xLower, opts2); tFs = toc(t0); fValFs = fCheck(rootFs); plotIteration(xx, ff, f, [xLower, rootFs], rootFs); itFs = 1; methodsData(iMethod).methodName = methodName; methodsData(iMethod).root = rootFs; methodsData(iMethod).fVal = fValFs; methodsData(iMethod).iterations = itFs; methodsData(iMethod).timeSec = tFs; catch tFs = toc(t0); methodsData(iMethod).methodName = methodName; methodsData(iMethod).root = NaN; methodsData(iMethod).fVal = NaN; methodsData(iMethod).iterations = NaN; methodsData(iMethod).timeSec = tFs; end else methodsData(iMethod).methodName = methodName; methodsData(iMethod).root = NaN; methodsData(iMethod).fVal = NaN; methodsData(iMethod).iterations = NaN; methodsData(iMethod).timeSec = NaN; end iMethod = iMethod + 1; %% 8) Method #7: vpaNewton (Symbolic) methodName = 'vpaNewton (Symbolic)'; if exist('vpa','file')==2 && license('test','Symbolic_Toolbox') figure('Name',methodName,'Color','w'); plot(xx, ff,'b-','LineWidth',2); hold on; grid on; xlabel('x'); ylabel('f(x)'); title(methodName); t0 = tic; [rootVpa, itVpa, xHistV] = vpaNewtonMethod(f, df, xLower, tol, maxIter); tVpa = toc(t0); rootV = double(rootVpa); fValVpa = fCheck(rootV); plotIteration(xx, ff, f, xHistV, rootV); methodsData(iMethod).methodName = methodName; methodsData(iMethod).root = rootV; methodsData(iMethod).fVal = fValVpa; methodsData(iMethod).iterations = itVpa; methodsData(iMethod).timeSec = tVpa; else methodsData(iMethod).methodName = methodName; methodsData(iMethod).root = NaN; methodsData(iMethod).fVal = NaN; methodsData(iMethod).iterations = NaN; methodsData(iMethod).timeSec = NaN; end iMethod = iMethod + 1; %% 9) Create table of results disp(' '); disp('------------------- RESULTS -------------------'); % We should have exactly 7 entries in methodsData now. methodNames = {methodsData.methodName}'; % 7x1 cell rootVals = [methodsData.root]'; % 7x1 double fVals = [methodsData.fVal]'; % 7x1 double iters = [methodsData.iterations]'; % 7x1 double times = [methodsData.timeSec]'; % 7x1 double % Format them as strings or keep numeric, then build table rootStrings = arrayfun(@(r) sprintf('%.6f', r), rootVals, 'UniformOutput',false); fValStrs = arrayfun(@(v) sprintf('%.8g', v), fVals, 'UniformOutput',false); iterStrs = arrayfun(@(i) sprintf('%d', i), iters, 'UniformOutput',false); timeStrs = arrayfun(@(t) sprintf('%.4g', t), times, 'UniformOutput',false); resultsTable = table(methodNames, rootStrings, fValStrs, iterStrs, timeStrs, ... 'VariableNames',{'Method','Root','f_at_root','Iterations','Time_sec'}); disp(resultsTable); end %==================================================================== % UTILITY #1: Evaluate f(x) if valid %==================================================================== function val = fCheck(x) % Return f(x) if x>0 and real, else NaN if x>0 && isreal(x) val = 10*log10(x) - x; else val = NaN; end end %==================================================================== % UTILITY #2: Plot iteration points & final root %==================================================================== function plotIteration(xx, ff, f, xHist, xFinal) plot(xx, ff,'b-','LineWidth',2,'DisplayName','f(x)'); hold on; grid on; yyHist = arrayfun(f, xHist); plot(xHist, yyHist,'ro-','LineWidth',1.5,'DisplayName','Iterations'); yFinal = f(xFinal); plot(xFinal, yFinal, 'ks','MarkerFaceColor','k','MarkerSize',8,... 'DisplayName','Root'); legend('Location','best'); hold off; end %==================================================================== % METHOD A: Bisection %==================================================================== function [root, iter, xHistory] = bisectionMethod(f, a, b, tol, maxIter) if f(a)*f(b) > 0 warning('Bisection: no sign change in [a,b].'); end xHistory = []; for i=1:maxIter c = (a+b)/2; xHistory(end+1) = c; %#ok fc = f(c); if abs(fc) < tol root = c; iter=i; return; end if f(a)*fc < 0 b = c; else a = c; end if abs(b-a) < tol root=c; iter=i; return; end end root=(a+b)/2; iter=maxIter; end %==================================================================== % METHOD B: Secant %==================================================================== function [root, iter, xHistory] = secantMethod(f, x0, x1, tol, maxIter) if x0<=0, x0=1e-6; end if x1<=x0, x1 = x0+1; end xHistory = [x0,x1]; f0 = f(x0); f1 = f(x1); for i=1:maxIter denom = (f1 - f0); if abs(denom)<1e-15 root=x1; iter=i; return; end x2 = x1 - f1*(x1 - x0)/denom; xHistory(end+1) = x2; %#ok if abs(x2 - x1) if abs(x1 - x0) 0 warning('FalsePosition: no sign change in [xl,xu].'); end xHistory = []; fl = f(xl); fu = f(xu); for i=1:maxIter xr = xu - fu*(xl - xu)/(fl - fu); xHistory(end+1) = xr; %#ok fr = f(xr); if abs(fr) end root=x0; iter=maxIter; end