SOLUTION OF EQUATIONS USING MATLAB
The present tutorial deals exclusively with numerical methods of solving algebraic and trigonometric equations, both single and several simultaneously, both linear and non-linear. Analytical solutions can be obtained using the methods described in the symbolic tutorial. When both symbolic and numerical methods work, sometimes one is simpler and sometimes the other.
In order to benefit from the following material one must copy and paste the indicated commands into MATLAB and execute them. These commands are shown in the format >> help plot, for example. The user should look up each command in MATLAB's help in order to understand the command and the different possible syntaxes that can be used with it.
Graphical solution of one non-linear equation or two non-linear equations in explicit form
When a single equation or a pair of equations is to be solved, it is good practice to first make a plot. In this way, it can be seen immediately what the approximate solution is.
In order to plot an equation using the "plot" command it must be in the explicit form, i.e. y = f(x). If the equation can be written only in implicit form, i.e. f(x,y) = 0, then the "ezplot" command must be used as described in the symbolic tutorial.
We consider first a single equation, i.e. f(x) = 0. As an example we will use sin2(x) e-x/2 – 1 = 0. Do the following steps:
>> x = -10:0.01:10;
gives x from -10 to +10 in increments of 0.01.
>> y = sin(x).^2.*exp(-x/2) -1;
Notice the dots before ^ and *. These dots are essential. Because x is a vector, we must use array operations rather than scalar operations (without the dot).
>> y0 = zeros(1,length(x));
Next let us consider a pair of explicit equations, y = sin2(x) e-x/2 – 1 and y = 10sin(x)/x. Proceeding as above:
>> clear; x = -10:0.01:10; y1 = 10*sin(x)./x; y2 = sin(x).^2.*exp(-x/2) - 1; plot(x,y1,x,y2);
Again we see that there are many solutions (intersections), both positive and negative. (Note that when MATLAB calculated sin(0)/0 it gave a warning, but nonetheless completed the other calculations and the plot.) As above, find the approximate solution at x near 3.5. (I get x = 3.49, y = 0.9796.) Write down your result to compare with those to be found later. Alternately, we can equate these two equations (since both = y), move everything to the left-hand side, and find the values of x when this is 0. For the value of x found, the corresponding value for y is determined by substituting this back into either one of the original equations.
>> clear; x = -10:0.01:10; yeq = 10*sin(x)./x - sin(x).^2.*exp(-x/2)+1;
>> yeq0 = zeros(1,length(x)); plot(x,yeq,x,yeq0);
Using the same method as above, I again got x = 3.49 and yeq = 0. Note that yeq is not y. To find the y we substitute x back into both of the simultaneous equations, which gives us two results. (Use your own value of x rather than my value of 3.49.)
>> clear; x = 3.49; y1 = 10*sin(x)/x, y2 = sin(x)^2*exp(-x/2) - 1
The closer these are y1 and y2, the better our solution.
Numerical solution of one or two non-linear equations in explicit form (y = f(x))
We use the "fzero" command, which finds the value of x for f(x) = 0. (See >> help fzero.) We will use the same examples as above. To display the result to 14 significant figures we first change the format to long. For the single equation, sin2(x) e-x/2 – 1 = 0:
>> clear; format long; x = fzero('sin(x)^2*exp(-x/2)-1', -1)
(Notice that no dot is required now before ^ or *, although using the dot causes no problem.) How does this result compare to the approximate result you obtained by the graphical method above?
Now we again consider the pair of equations, y = sin2(x) e-x/2 – 1 and y = 10sin(x)/x:
>> clear; x = fzero('10*sin(x)/x - sin(x)^2*exp(-x/2) + 1', 3.5)
>> y1 = 10*sin(x)/x, y2 = sin(x)^2*exp(-x/2) - 1
As before, we see that the numerical method is more accurate. So why bother with the graphical method at all? The reason is, to see about where you need to have "fzero" seek a solution!
Finding a particular solution when there are infinitely many
There are equations with an infinite number of solutions, for example sin(x) = 1/2.
It is helpful to see some of these solutions by plotting y = sin(x) - 1/2 and then observing where y has values of zero:
>> clear, syms x; eq='y - sin(x) + 1/2'; ezplot(eq,[-6,6,-2,1])
>> hold on; eq0='0'; ezplot(eq0); hold off
The "hold on" command tells MATLAB to write the subsequent plots on top of the first one, rather than replacing it. Plotting 0 puts a horizontal line on the graph. Intersections of the sine wave with this line represent solutions of the first equation. Approximate values can be obtained in the Figure window by clicking on Tools / Zoom In, clicking on a desired intersection to enlarge that area of the graph, clicking on Tools / Data Cursor, and then clicking on that intersection again to give the approximate values of x and y there. Try this on the first two intersections to the right of the origin (x > 0). Write down your results to compare with those found below.
Here are the steps to find a solution numerically:
function out = func2(x)
out = sin(x) - 1/2;
>> clear, x = - 4:0.01:4; plot(x,func2(x))
>> fzero('func2',3) or >> fzero(@func2,3) or >> fzero('func2', [2,3])
MATLAB finds the value of x nearest 3 that gives func2 = 0. Notice (“whos” or Workspace) that the result is double, so no conversions are necessary for subsequent numeric calculations.
Numerical solution of two or more equations in implicit forms (e.g., f1(x,y)=0, f2(x,y)=0)
We use the "fminsearch" command. Consider as an example the following two implicit equations (already in MATLAB format):
0.5 = (200+3*x+4*y)^2 / ((20+2*x+3*y)^2 * x)
10 = (20+2*x+3*y)*y / x
Use the following steps:
% function to calculate the values of x and y
% satisfying the two equibiliubrium relations
% for a methanol reactor, CO + 2H2 = CH30H
% CO2 + 3H2 = CH3OH + H2O (coal to chems case study)
% x is the molar flow rate of CO, y of CO2 (kmol/h)
% After saving: >> fminsearch('meth_reac',[25,3]); x=ans(1), y=ans(2)
% Any positive values work as initial guesses in this example
eq1 = 0.5 -(200+3*x+4*y)^2/(20+2*x+3*y)^2/x;
eq2 = 10 -(20+2*x+3*y)*y/x;
3. In the Command window, type "fminsearch('function name',[guessed values]). For our example,
>> fminsearch('meth_reac',[25,3]); x=ans(1), y=ans(2)
You can check your result by using the symbolic "solve" command, as follows for our example:
>> clear, syms x y
>> [x y]=solve(eq1,eq2,x,y)
Note that several solutions are produced, from which you have to select the physically reasonable one. If this is the first solution, then to obtain double-precision numerical variables we use:
>> x=double(x(1)), y=double(y(1))
(The equations in this example are based on the equilibrium relationships for a chemical reactor making methanol from CO, CO2 and H2 as part of an existing plant that produces a variety chemicals starting with the gasification of coal, i.e. coal reacting with water at high temperature.)
Numerical solution of simultaneous linear equations using the matrix backslash (\) method
MATLAB uses Gaussian elimination to solve systems of linear equations via the backslash matrix command (\). As an example, consider the following three equations:
The procedure is as follows:
2x - 3y + 4z = 5
x + y + 4z = 10
3x + 4y - 2z = 0
>> clear, C = [2,-3,4; 1,1,4; 3,4,-2], B = [5; 10; 0]
>> A = C\B, x = A(1), y = A(2), z = A(3)
Check your results using: >> C*A - B
(Note that these are matrix operations; you MUST NOT use a dot (.) before * or \)
Compare this result with that obtained using MATLAB symbolic mathematics for the same system of equations.
Another application of the backslash matrix operator is to fit data to a non-linear equation. We use an example with the following y versus data:
>> clear; t = [0, .3, .8, 1.1, 1.6, 2.3]'; y = [0.5, 0.82, 1.14, 1.25, 1.35, 1.40]'; plot(t,y,'o')
Notice that this does not give a straight line plot. Instead we will try a quadratic fit (without using the Basic Fitting tool on the graph window or the "polyfit" command). That is, we want to find the values of the coefficients (a0, a1, a2) that best fit these data to the form y = a0 + a1t + a2t2. We have 6 values each of t and y, from which we generate 6 equations in the format a0 + ta1 + t2a2 = y matching that required for solution of simultaneous linear equations:
a0 = 0.5
a0 + 0.3a1 +0.32a2 = 0.82
a0 + 0.8a1 +0.82a2 = 1.14
a0 + 1.1a1 +1.12a2 = 1.25
a0 + 1.6a1 +1.62a2 = 1.35
a0 + 2.3a1 +2.32a2 = 1.40
Thus we have six simultaneous equations, but only 3 unknowns. We cannot expect to get values for a0, a1 and a2 that would exactly satisfy all 6 equations. The backslash operator gives the least squares values instead. We must put these 6 equations in the form C*A=B, as follows (assuming t and y have already been entered, as above):
>> C(:,1) = ones(length(t),1)
>> B = y
>> A = C\B
>> a0 = A(1), a1 = A(2), a2 = A(3)
>> tfit=0:0.1:2.3; yfit = a0 + a1*tfit + a2*tfit.^2;
This is a general method and not confined to fitting data to polynomials. See fitting data to non-polynomial equations.