%% Example of variety of ways of finding roots of a nonlinear equation.

% NOTE: there are several "pause" commands below - program suspends
% until a key is pressed - this is to allow you to view a stage in the
% program before proceeding to the next stage.

% This example solves Van der Waal's equation for molar volume (v) given
% pressure (p) and temperature (T). This example uses the inputs of
% Himmelblau & Riggs, "Basic Principles and Calculations in Chemical Engineering,"
% 7th ed., example 15.12, p. 465
% with the CHANGE that the number of moles is fixed at n = 1 lb-mol here.
%
% Van der Waals equation of state for fluids can be arranged several ways
% Two of the three state variables are easy to obtain:
%  (p + a/v^2)*(v - b) = R*T  >> given v and p, solve for T
%   p = -a/v^2 + R*T/(v - b) >> given v and T, solve for p
% Solving for v given p and T is more complex and is the problem we consider
%   v^3 - (b+R*T/p)*v^2 + a*v/p - a*b/p = 0
% A problem example is, given temperature and pressure and vessel volume,
% find the number of moles of a fluid that can be contained in the vessel.
%
% This nonlinear equation happens to be a 3rd order polynomial.
% The standard Matlab function roots() used below works only with polynomials.
% All the other methods illustrated below work with most any type of
% nonlinear equation, not just polynomials.

fprintf('---------------------------\n') % run separator
clear all

% Given: conditions in vessel
p = 679.7;  % pressure (psia)
T = 683; % (deg R)

% Given: Van der Waal's equation coefficients
a = 3.49E4;
b = 1.45;
R = 10.73; % ideal gas constant (psia-ft3)/(lbmol-R)

% vector of molar volume values to use in computations for plot
v = [0:0.1:10];  

% note need array operator (.^) with v since v is an array here
f = v.^3 - (b+R*T/p)*v.^2 + a*v/p - a*b/p;

z = 0*v;

plot(v,z,'g',v,f,'b')

xlabel('v, molar volume of propane, ft3')
ylabel('f(v) = v.^3 - (b+R*T/p)*v.^2 + a*v/p - a*b/p')
title('see Himmelblau & Riggs, 7th ed., Ex 15.12, p. 465')

fprintf('\nThis shows several different techniques to find the root of a nonlinear equation. \n')
fprintf('The equation we want to find the root of is f(v) = v.^3 - n*(b+R*T/p)*v.^2 + n^2*a*v/p - n^3*a*b/p \n\n')
fprintf('******* GO TO THE PLOT WINDOW AND SEE IF YOU CAN DETERMINE THE ROOT FROM THE PLOT ******* \n\n')

fprintf('CLICK IN COMMAND WINDOW & PRESS ANY KEY TO CONTINUE \n\n')
pause

%% use roots() function

fprintf('******* Now use roots() function of MATLAB ******* \n\n')
fprintf('Since this equation is a polynomial, we can use the roots() function of MATLAB \n')
fprintf('c vector holds the 3rd-order polynomials coefficient values [1, -(b+R*T/p), a/p, -a*b/p] \n')

c = [1, -(b+R*T/p), a/p, -a*b/p]

fprintf('We want the positive real root \n')
theRoots = roots(c)

fprintf('roots() is done \n') 
fprintf('CLICK IN COMMAND WINDOW & PRESS ANY KEY TO CONTINUE \n\n')
pause

%% use Newton's method

fprintf('******* Now doing Newtons method ******* \n\n')

endTol = 1e-4; % see line below for getting this from user
% endTol = input('Input the tolerance, fractional change between current and previous est., e.g., 1e-4 (try others) \n');

v = R*T/p  % initial guess using v estimated from ideal gas law

% note v now becomes a scalar

thisTol = 100;  % fractional change in v between iterations
                % set initial value high so goes into repeat at least once
				
% use a counter so we only display v value a few times
c = 0;

while thisTol > endTol

    oldv = v; % need to save old v value for thisTol calculation below
	
    % slope = df/dv = ( f(old v) - 0 ) / (old v - new v)
    % solve for new v = old v - f(old v)/slope
     
    f = v^3 - (b+R*T/p)*v^2 + a*v/p - a*b/p;
    slope = 3*v^2 - 2*(b+R*T/p)*v + a/p;
    v = v - f/slope;
    
    % here we used an analytical expression for df/dv
    % For more complex equations we can estimate df/dv numerically at each
    % step by computing f at two different and relatively closely spaced values of v
	
    thisTol = abs((v - oldv)/v);
    
    c = c+1;
    if c < 6
        v % echo value to screen
    end
	
end

fprintf('Newtons method estimate, v = %6.4f \n\n',v)
fprintf('CLICK IN COMMAND WINDOW & PRESS ANY KEY TO CONTINUE \n\n')
pause

%% using fzero() function

fprintf('******* Now use fzero() function of MATLAB ******* \n\n')
fprintf('To specify the function in one line instead of external function file, we need to \n')
fprintf('use the variable name x instead of v and we have to supply the coefficient values as constants. \n')
fprintf('By defining and external function file, we can have the computer calculate coefficients. \n\n')

v = fzero('x^3 - 12.232*x^2 + 51.246*x - 74.452', [0 10], 1e-4, 1)

fprintf('fzero() is done \n\n')
fprintf('CLICK IN COMMAND WINDOW & PRESS ANY KEY TO CONTINUE \n\n')
pause

%% using fminbnd() function

fprintf('******* Now use fminbnd() function of MATLAB ******* \n\n')

fprintf('Since fminbnd() finds x that minimizes f(x), and since our original f can go \n')
fprintf('negative, we need to find the minimum of the absolute value, or of the square of f \n\n')

v = fminbnd('abs(x^3 - 12.232*x^2 + 51.246*x - 74.452)', 0, 10, [1, 1e-4])

fprintf('fminbnd() is done\n')
fprintf('For functions of two or more adjustable variables, use fminsearch() \n\n')
fprintf('CLICK IN COMMAND WINDOW & PRESS ANY KEY TO CONTINUE \n\n')
pause

%% using Interval Halving method

fprintf('******* Now use the interval halving (bisection) method ******* \n\n')

% specify interval over which f(v) changes sign
Vmin = 2;
Vmax = 8;

% use a counter so we only display v value a few times
c = 0;

fVmid = 1; % arbitrary number so do repeat at least once
% repeat is not optimized for speed here
while abs(fVmid) > 0.0001
    
    fVmin = Vmin^3 - (b+R*T/p)*Vmin^2 + a*Vmin/p - a*b/p;
    fVmax = Vmax^3 - (b+R*T/p)*Vmax^2 + a*Vmax/p - a*b/p;
    
    Vmid = (Vmin + Vmax)/2;
    fVmid = Vmid^3 - (b+R*T/p)*Vmid^2 + a*Vmid/p - a*b/p;
    
    c = c+1;
    if c < 6
        Vmid % echo value to screen
    end
    
    if sign(fVmid) == sign(fVmax)
        Vmax = Vmid;
    else
        Vmin = Vmid;
    end
	
end

fprintf('Bisection Method estimate, v = %6.4f \n\n',Vmid)
fprintf('CLICK IN COMMAND WINDOW & PRESS ANY KEY TO CONTINUE \n\n')
pause

%% using Successive Approximation method

fprintf('******* Now use the successive approximation method ******* \n\n')

v = [0:0.1:10];  % specify array of v values to make a plot

% now v is an array so need to use array (dot) operators .* and and .^

% now the function f(v) is arranged so f(v) = v, rather than = 0
fv = -(p/a)*(v.^3 - (b+R*T/p).*v.^2 - a*b/p);

plot(v,fv)
hold on

plot(v,v)

title('Plot for Successive Approximation method, Himmelblau & Riggs, 7th ed., Ex 15.12, p. 465')
xlabel('v, molar volume of propane, ft3')
ylabel('fv = -(p/a)*(v.^3 - (b+R*T/p).*v.^2 - a*b/p)')
hold off

fprintf('******* GO TO THE PLOT WINDOW AND SEE THE PLOT FOR SUCCESSIVE APPROXIMATION ******* \n\n')

v = R*T/p; % initial guess using ideal gas law, v is a scalar now

thisTol = 100;  % fractional change in v between iterations
                % set initial value high so goes into repeat at least once

while thisTol > 1e-4 % ending tolerance on changes in v values
	
	vlast = v;
	v = -(p/a)*(v^3 - (b+R*T/p)*v^2 - a*b/p);
	thisTol = abs((v - vlast)/v);

end

v % echo final v to command window

fprintf('only limited number of iterations performed \n\n')
fprintf('CLICK IN COMMAND WINDOW & PRESS ANY KEY TO CONTINUE \n\n')
pause

%% use SYMBOLIC TOOLBOX

fprintf('******* Now use the SYMBOLIC TOOLBOX ******* \n\n')

% ALSO SEE LINK TO WOLFRAM ALPHA NEAR TOP OF MY MATH TOOLS PAGE
% WOLFRAM ALPHA CAN ALSO DO SYMBOLIC MATH FOR YOU

fprintf('First solve with all as symbols - see comments in program \n')
fprintf('   want the first line of answer (scroll right), which is the real root \n')
fprintf('   copy it, paste into command window, it evaluates to 4.4844 \n')

% SOLVE() WILL SOLVE FOR x SO SAFEST TO USE x
% BUT IF x DOES NOT APPEAR, AS IT DOES NOT HERE, SOLVE() WILL
% ATTEMPT TO SOLVE FOR VARIABLE NEAREST IN ALPHABET, WHICH IS v HERE

solve('v^3 - (b+R*T/p)*v^2 + a*v/p - a*b/p')
 
% the first line of ans is a real expression of constants
% the line below was copied from command window as ans from line above

((R*T + b*p)^3/(27*p^3) + (((R*T + b*p)^3/(27*p^3) - (a*(R*T + b*p))/(6*p^2)...
    + (a*b)/(2*p))^2 - ((R*T + b*p)^2/(9*p^2) - a/(3*p))^3)^(1/2) - ...
    (a*(R*T + b*p))/(6*p^2) + (a*b)/(2*p))^(1/3) + (R*T + b*p)/(3*p) ...
    + ((R*T + b*p)^2/(9*p^2) - a/(3*p))/((R*T + b*p)^3/(27*p^3) + ...
    (((R*T + b*p)^3/(27*p^3) - (a*(R*T + b*p))/(6*p^2) + (a*b)/(2*p))^2 ...
    - ((R*T + b*p)^2/(9*p^2) - a/(3*p))^3)^(1/2) - (a*(R*T + b*p))/(6*p^2)...
    + (a*b)/(2*p))^(1/3)

% copy this and paste into command window and evaluate and get
% 4.4844, which is the value of the real root returned above by roots()

fprintf('Next, solve with coefficients as constants - see comments in program \n')
fprintf('   want the first line of answer, which is the real root \n')

% evaluate coefficients as constants
c1 = (b+R*T/p); % = 12.2321 for original parameter values
c2 = a/p; % = 51.3462
c3 = a*b/p; % = 74.4520

solve('v^3 - 12.2321*v^2 + 51.3462*v - 74.4520')

fprintf('For symbolic math, also see link to Wolfram Alpha examples near top of my Math Tools page. \n\n')

fprintf('Note how different methods gave different answers. \n')
fprintf('The same method can give different answers with different initial guesses and tolerances. \n\n')
fprintf('PROGRAM HAS ENDED  \n\n')