Download Numerical Methods Lecture 3: Solving Nonlinear Equations and Finding Roots and more Study notes Civil Engineering in PDF only on Docsity!
Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods
Lecture covers two things:
- Solving systems of linear equations symbolically
- Using Mathcad to solve systems of nonlinear equations
- Investigating algorithms to find roots of equations
Solving Systems of Linear Equations Symbolically
Let’s take a look at a very powerful tool in Mathcad that started a revolution in computational analysis. It started in the late 1980’s when I was an undergraduate. A company called Wolfram created a computer program called Mathematica. This was the first computer code that could solve algebraic and calculus equations symbolically. That is, if I had an equation that said x*y = z, Mathematica could tell me that y = z / x, without ever needing me to assign numbers to x, y, or z. It also was able to solve integrals, dif- ferential equations, and derivatives symbolically. This was an incredible advance, and opened the doors to a whole new world of programming, numerical methods, pure mathematics, engineering, and science.
Since then, a competing code called Maple was developed and sold itself to other software companies to include in their programs.
The end result: Mathcad uses Maple as a solving engine in the background (you don’t see it) to solve problems symbolically. Here we will look at a brief example of how to use this capability in the context of solving a system of linear equations.
Example : The structural system below is something you will see in CES 3102 or CES 4141. r1 and r2 are labels that indicate how the ends of the beam are allowed to move. Q and W represent the external loads (a couple and a distributed load, respectively), and material properties are given as E and I. L is the length of the beam. The goal is to solve for the amount or rotation at r1, and the deflection at r2 that occurs for thegiven loads. This would help us to solve for internal stresses, allowing us to design the beam to sur- vive these internal forces.
Solution : The way we learn to solve this problem in CES 4141 is using a Matrix-based solution proce- dure. The generic form of the solution is K * r = R K is a 2 x 2 ‘stiffness’ matrix that contains information about the structures shape, boundary conditions and material properties. THe information needed is L, E, and I
r is a 2 x 1 vector that contains the unknown rotation and displacement quantities sought.
R is a 2 x 1 vector that contains only information about the external loads (W and Q for this problem).
r1 r
W
Q
L
W = 3 K/FT
Q = 2 K*FT
L = 15 FT
E = 4000 ksi I = 1400 in^
Solution continued : So we will have K as a known 2 x 2 stiffness matrix We will have R as a known 2 x 1 load vector We will solve for the unknown displacement vector r
We will not go into any detail on HOW we fill in the values for K and R, that’s for another class. Here we will just consider how to solve the ssytem of matrix equations symbolically.
Define Load Vector R as a Define Stiffness K as a function of E, I, L function of Q, W, L
K E I( , ,L)
4 E⋅ ⋅I
L
6 E⋅ ⋅I
L
2
6 E⋅ ⋅I
L
2
12 E⋅ ⋅I
L 3
:= R Q W(^ , ,L)
Q
W L
2 ⋅ 30
⋅ W⋅L
The variables in the parenthesis (E,I,L) and (Q,W,L) is telling mathcad that the variables K and R, respectively, will be a function of those variables.
known information
Show me the inverse of the stiffness matrix Not a necessary step, but pretty cool see the matrix version of 1/K Note the arrow --> is the way to initiate a symbolic operation. Here we are asking what is the inverse of K? The arrow comes from the view => toolbars => symbolic pad
K E I( , ,L) 1 −
1 ( E I⋅)
⋅L
− 1 ( 2 E I⋅ ⋅)
⋅L^2
− 1 ( 2 E I⋅ ⋅)
L ⋅^2
1 ( 3 E I⋅ ⋅)
⋅L^3
→
Calculate the displacement vector as K -1^ * R Note that we DO need to keep writing the names of the independent variables next to K, r, and R to tell Mathcad what symbols to look for
r E I( , , L, Q,W) K E I( , ,L) 1 − := ⋅R Q W( , ,L)
Display the resulting displacement vector calculated in terms of W, Q, L, E, I
Note the arrow displays a r E I( , , L, Q,W) symbolic result
1 (E I ⋅)
⋅ (^) L Q 1 30
⋅ (^) WL
− ⋅^2
⋅ 7 ( 40 E I⋅ ⋅)
L ⋅ 3 − ⋅W
− 1 ( 2 E I ⋅ ⋅)
L ⋅ 2 Q
1 30
⋅W (^) L
− ⋅^2
⋅ 7 (60 E I ⋅ ⋅)
L ⋅ 4
→
Symbolic Solution
Symbolic Inverse
Symbolic Calculation
Solving Systems of Nonlinear Equations
We won’t go into the algorithms themselves here. We will just focus on how to use Mathcad to solve the problem.
Use Mathcad help and use the keywords ‘nonlinear equations’ to get some information. Follow the link to the ‘Find’ function in Mathcad. ‘Find’ is the workhorse that gets the solution. All you need to do is give it the correct description of the problem. This includes the equations to solve, any possible constraints to the solutions you are interested in, and a place to start looking for a solution (i.e., initial guess for the solution). Here we will go straight into an example to show how it’s done.
f3 x( ) (^) := 1 x .15 x
6 4 2 0 2 4 6
4
2
0
2
4
6
f1 i( ) f2 i( ) f3 i( )
i
ORIGIN ≡ 1
Using the Find function to solve systems of nonlinear equations This example will solve for the intersecting values of the followin system of 2 equations
x^2 + (y − 2 ) 2 8
x + .15 x ⋅^2 + y 1
Let's first take a look at the solution visually by plotting We will create some functions that allow us to plot the two equations above
f1 x( ) := 8 − x^2 + 2
i :=− 5 , −4.99.. 4 f2 x( )^ :=−^8 −x^2 +^2 f3 x( ) := 1 − x−.15 x ⋅^2 f3( x ) (^) := 1 x .15 x
6 4 2 0 2 4 6
4
2
0
2
4
6
f1 i( ) f2 i( ) f3 i( )
i
We can see that there are two solutions to the two equations above. We'll try to find both solutions by using the 'find' function
Use the help desk with the keyword 'find function'
Step one: start by guessing the solution to the set of two equations we are trying to solve. Let's guess that x=4 and y= -
x :=^4 y :=− 1
Next we create what is called a 'solve block'. This defines the equations to be solved.
- type the word 'Given'
- type below 'Given' all the equations to be solved
Note that we use the 'boolean' tab to get the equal sign, (^) not the keyboard stroke.
Given
x^2 + (y − 2 ) 2 8
x + .15 x ⋅^2 + y 1
Finally, we use the 'Find' function to solve by sending in the two guesses we made above for x and y
Find x y( , )
−0.
=
'Find' starts looking at the initial guess values, and iteratively updates the values of x and y until a pair is found that solves the equations defined in the solve block.
The results give the x and y values corresponding to one of the solutions. find where x=1.278 and y= -0.523 in the graph above 'Find' found this solution because the initial guesses were closer to this solution.
Solving for Roots of Equations
What?
- Locating where some function crosses the x-axis Find x such that
Why?
- We can then find where crosses any value
Find x such that or Find x such that
When?
- Numerical methods in general are:
- less accurate
- slower than using an analytical (exact) solution.
- Use numerical methods when analytical answers not available
Example:
Given: , Find x such that y = 0
analytical solution: ==> no need for
Given: , Find x such that y = 0
No idea? Good, a prime candidate for a numerical solution
We will consider several methods, each does the following in some way:
- Start with some initial guess (or range in which to look)
- Use some information about the function to make another, closer guess
- Stop when our current guess is ‘close enough’ to a solution
The user will decide what is ‘close enough’ (how big the error can be before we can stop looking)
Method 1: INCREMENTAL
SEARCH METHOD
x
y
x 1 first guess
error (y distance from 0)
b a
search in the range [a,b] for a root
x
y
x 2 next guess
new error (y distance from 0)
b a
Start with an initial range that contains the root, and subdivided it into several smaller sub-ranges.
Look inside each sub-range one by one for the root. When the sub-range containing the root is identi- fied, choose either end of the range as the guess.
Evaluate the error in your guess. If error is too big, subdivide your new range into smaller sub-ranges.
Loop back to step 1) and continue the loop until the error in step 3) is small enough
- Out of the 4 subranges, we pick the one in which the function crosses the x-axis as the new range to subdivide even smaller.
- Visually we can see this, but how does the algorithm know how to identify the smaller range?
- Remember we have the function, so we can evaluate it at any x location we choose.
Take a moment, work out a way to identify the new smaller subrange.
a 1
b 1
a 2
b
start looking in range [a 1 , b 1 ] segment into 4 pieces (we picked 4 out of thin air)
look in each subrange until [a 2 , b2] is found
Now divide [a 2 , b2] into 4 segments and look again
evaluate the error err = minimum ( | f(a 2 ) |, | f(b 2 ) | )
subrange subrange
subrange
subrange
a 2
b a 3
b 3
look in each subrange until [a 3 , b3] is found evaluate the error err = minimum ( | f(a 3 ) |, | f(b 3 ) | )
subrange subrange
subrange
subrange
Now divide [a 3 , b3] into 4 segments and look again ...
Method 2: BISECTION METHOD
Start with an initial range [a, b], and cut the range in half, assume this half point m is the root
Evaluate the error err = | f(m) |.
If error is too big, decide if the root is to the left or right of [m].
- If root is to the left of m, reassign a new range a = a, b = m
- Else If root is to the right of m, reassign a new range a = m, b = b
- Go back to step 1), stop bisecting when the error at m is ‘small enough’.
Three Bisection iterations
- Eventually we zero in on the root
- How do we decide where the root is relative to m for the next [a , b] ??
Same as the incremental search:
prod = f(a) * f(m)
if prod < 0, root is to the left
otherwise root is to the right
a
m b
err > tol - yes?
evaluate the error err = | f(m) |
a
err = |f(m)|
err = |f(m)|
new m
make new range: root is to left of m ===> a = a, b = m root is to right of m ==> a = m, b = b divide new [a , b] in half to get new m evaluate the error err = | f(m) | (^) ....etc.
new b
a
m1 b
f(m1)
iteration 1
b m
f(m2)
a
iteration 2
a
m3 b
f(m3)
iteration 3
b stays, a moves over from left
error at m1 too big, so guess again
error at m2 too big, so guess again
- This slow method does not use much function information. We can do better.
Method 3: FALSE POSITION
- Relative height of function at end points used to make better guesses
- Define initial range [a b] (possibly the result of a single pass of the incremental search method).
first guess c - intersection of a line from f(a) to f(b) with x-axis
(from similar triangles)
- Evaluate the error err = | f(c) |
If error is too big, decide if the root is to the left or right of [c]. for f(a) * f(c) < 0, reassign a new range a = a, b = c for f(a) * f(c) > 0, reassign a new range a = c, b = b
new guess
- Stop iterating when the error is small enough
b
a
f(a)
f(b)
it seems more likely that the root would be closer to a than b, and not in the middle
Given the heights of f(a) vs. f(b)
Let’s use this concept in an algorithm
this is from the definition of slope
xk x
tangent to f(xk)
xk+
solve for x (^) k+
f(xk)
Let’s continue with a few more iterations
Another example, 2 iterations in one figure
x 1 x first guess
error (y distance from 0)
tangent to f(x 1 ) second guess x 2
f(x 2 )
x 1 x
tangent to f(x 2 )
second guess x (^2)
f(x 2 )
x (^3)
x 1 x
x (^2)
f(x 3 )
x 3
x
x
f(x1)
slope = f’(x1)
x
f(x2)
slope = f’(x2)
x
