*W.R. Wilcox, Clarkson
University*

*Last
revised May 13, 2010*

**SOLUTION OF****
EQUATIONS USING MATLAB**

**See also: **

·
**Tutorial on
symbolic mathematics using MATLAB**

·
**Relationships
between mathematical, MATLAB and Excel expressions**

·
**MATLAB tutorials
on analysis and plotting of data**

·
**Common conversion
factors for units**

·
**Introductory
MATLAB tutorials: See MATLAB help, MATLAB, getting started. The following sites also have useful
tutorials: Mathworks, Florida, Louisiana,
New
Hampshire, Carnegie
Mellon**

**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 sin ^{2}(x)
e^{-x/2} – 1 = 0. Do the
following steps:**

**While optional, it's a good idea to clear the Workspace memory in order to avoid the chance of MATLAB using a previous value for a variable:**

** >>
clear**

**Guess the range of x values over which you expect a solution and generate values for the x vector over this range. Select an increment small enough to yield a smooth curve. For example:**

** >>
x = -10:0.01:10;**

** gives x from -10 to +10 in increments of 0.01.**

**Generate the y = f(x) vector corresponding to these x values:**

** >>
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).**

**Also create a line for y = 0. The intersection of this line with the curve above gives the solution (or solutions, as in this example). To do this we create a vector with the same number of values as x, with each value being 0.**

** >>
y0 = zeros(1,length(x));**

**Make the plot:**

** >> plot(x,y,x,y0);**

**Locate the solution(s). If the two lines do not intersect, repeat with another range for x. In our example, the lines intersect more than once, and we surmise, in fact, that there are infinitely many solutions for negative x. Let us select the solution at about x = -1. Enlarge this area of the intersection in the Figure window by Tools / Zoom In, and then clicking near the intersection once or more. Then go to Tools / Data Cursor. Click as close as you can to the intersection. Write the result down to compare with that obtained by using the "fzero" command below.**

**Next let us consider a pair of explicit
equations, y = sin ^{2}(x) e^{-x/2} – 1 and y =
10sin(x)/x. **

** ****>> 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, sin ^{2}(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 = sin ^{2}(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:**

**Convert the equation to a function consisting of everything moved to the left-hand side, e.g. func2(x) = sin(x)-1/2.****Create this function in the MATLAB editor, save to the Current Directory, and make certain the current directory is in the path (File / Set Path). For example:**

**function****
out = func2(x)**

**out**** = sin(x) - 1/2;**

**end**

**Plot this function over the range of interest to see where the solutions are, e.g.:**

** >> clear, x = -
4:0.01:4; plot(x,func2(x))**

**Use one of the following formats to find the solution:**

** >>
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:**

**Move everything to the left-hand side in all equations.****Using MATLAB's edit window, create a function that calculates the deviations from 0 when arbitrary values are used for the variables. The final output of the function is the sum of the squares of these deviations. Following is the code for our example, and would be saved as meth_reac.m in MATLAB's Current Directory, which must also be in the Path (File/Set Path).**

**function**** out=meth_reac(guesses)**

**% 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**

**x=guesses(1); y=guesses(2);**

**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;**

**out=****eq1^2+eq2^2;**

**end**

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

**>> 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'**

**>> [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, CO _{2} and H_{2} 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:**

** 2x-3y+4z=5**

** y+4z+x=10**

** -2z+3x+4y=0**

**The procedure is as follows:**

**Each equation must first have each variable in the same order on the left-hand side, with the constants on the right-hand side. This should be done using pencil and paper. Students commonly assume they can do this while entering the information in MATLAB, and frequently make mistakes doing so. For our example, we write:**

**2x - 3y + 4z = 5**

**x + y + 4z = 10**

**3x + 4y - 2z = 0**

**In MATLAB, create matrices of coefficients for the variables and constants on the right-hand sides. Thus for our example:**

** >> clear, C = [2,-3,4; 1,1,4; 3,4,-2], B
= [5; 10; 0]**

**Solve for x, y and z using:**

** >> 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 (a_{0},
a_{1}, a_{2}) that_{ }best fit these data to the form y
= a_{0} + a_{1}t + a_{2}t^{2}. We have 6 values each of t and y, from which
we generate 6 equations in the format a_{0} + ta_{1} + t^{2}a_{2}
= y matching that required for solution of simultaneous linear equations:**

**a _{0} = 0.5**

**a _{0} + 0.3a_{1 }+0.3^{2}a_{2} = 0.82 **

**a _{0} + 0.8a_{1 }+0.8^{2}a_{2} = 1.14 **

**a _{0 }+ 1.1a_{1 }+1.1^{2}a_{2 }=
1.25 **

**a _{0 }+ 1.6a_{1 }+1.6^{2}a_{2 }= 1.35 **

**a _{0} + 2.3a_{1 }+2.3^{2}a_{2} = 1.40**

**Thus we have six simultaneous equations, but only 3
unknowns. We cannot expect to get values
for a _{0}, a_{1} and a_{2} 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)**

**>> C(:,2)=t**

**>> C(:,3)=t.^2**

**>> 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;**

**>> plot(t,y,'o',tfit,yfit);**

**This is a general method and not confined to fitting data
to polynomials. See ****fitting data to non-polynomial equations.**** **