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

1. 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

2. 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.

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

4. 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));

5. Make the plot:

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

6. 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²(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, sin²(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²(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:

1. Convert the equation to a function consisting of everything moved to the left-hand side, e.g. func2(x) = sin(x)-1/2.
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

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

4. 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:

1. Move everything to the left-hand side in all equations.
2. 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₂ and H₂ 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:

1. 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

2. 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]

3. 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₀, a₁, a₂) that best fit these data to the form y = a₀ + a₁t + a₂t². We have 6 values each of t and y, from which we generate 6 equations in the format a₀ + ta₁ + t²a₂ = y matching that required for solution of simultaneous linear equations:

a₀ = 0.5

a₀ + 0.3a₁ +0.3²a₂ = 0.82

a₀ + 0.8a₁ +0.8²a₂ = 1.14

a₀ + 1.1a₁ +1.1²a₂ = 1.25

a₀ + 1.6a₁ +1.6²a₂ = 1.35

a₀ + 2.3a₁ +2.3²a₂ = 1.40

Thus we have six simultaneous equations, but only 3 unknowns. We cannot expect to get values for a₀, a₁ and a₂ 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.