Spring 2024-2025
Numerical Methods for Science and Engineering
Lecture Note 4
Nonlinear System of Equations
4.1 Newton-Raphson Method
Recall that the Newton-Raphson formula to find the root of the equation 𝑓(𝑥) = 0 can be
written as
𝑥𝑛+1 = 𝑥𝑛 − [𝑓 ′ (𝑥𝑛 )]−1 𝑓(𝑥𝑛 ).
For a system of nonlinear equations, the Newton-Raphson formula in matrix form can be
expressed as
𝑋𝑛+1 = 𝑋𝑛 − [𝐽(𝑋𝑛 )]−1 𝐹(𝑋𝑛 ).
where X is a variables matrix, F(X) is a function matrix and J(X) is a Jacobian matrix.
The variable, function and Jacobian matrices for a system of two variables can be written
as
𝜕𝑓1 𝜕𝑓1
𝑥 𝑓1 (𝑥, 𝑦) 𝜕𝑥 𝜕𝑦
𝑋 = ( ), 𝐹=( ), 𝐽(𝑥, 𝑦) = (𝜕𝑓 ).
𝑦 𝑓2 (𝑥, 𝑦) 2 𝜕𝑓2
𝜕𝑥 𝜕𝑦
This method can be extended for equations with more variables.
𝑥2 + 𝑦 − 1 = 0
Example 4.1 Consider the nonlinear system .
𝑦 3 − 2𝑥 + 4 = 0
(a) Find the Jacobian matrix 𝐽(𝑥, 𝑦) for the above system.
(b) Evaluate the inverse of the Jacobian matrix at (1.4, −1).
(c) Write down the iterative formula for the above system based on the Newton-
Raphson method.
(d) Estimate the root to 2 d.p. using the using the above iterative formula once
starting with (𝑥0 , 𝑦0 ) = (1.4, −1).
(e) Write down MATLAB commands to execute the iteration four times.
(f) Use MATLAB function “fsolve(fun,x0)” to find the above root.
Solution:
(a) The function matrix and Jacobian matrix are
𝑥2 + 𝑦 − 1 2𝑥 1
𝐹(𝑥, 𝑦) = ( ), 𝐽(𝑥, 𝑦) = ( )
𝑦 3 − 2𝑥 + 4 −2 3𝑦 2
1
Spring 2024-2025
(b)
𝑎 𝑏
Note: The inverse of a 2 × 2 matrix 𝐴 = ( ) can be calculated by
𝑐 𝑑
1 𝑑 −𝑏
𝐴−1 = ( )
|𝐴| −𝑐 𝑎
At the point (1.4, −1). we have
2.8 1
𝐽(1.4, −1) = ( )
−2 3
and
[𝐽(1.4, −1)]−1 = (0.2885 −0.0962
)
0.1923 0.2692
(c) The iterative formula is
𝑋𝑛+1 = 𝑋𝑛 − [𝐽(𝑋𝑛 )]−1 𝐹(𝑋𝑛 ).
where
𝑥𝑛
𝑋𝑛 = (𝑦 ).
𝑛
(d) Using 𝑛 = 0, we have
𝑋1 = 𝑋0 − [𝐽(𝑋0 )]−1 𝐹(𝑋0 )
where
1.4 −0.04
𝑋0 = ( ), 𝐹(𝑋0 ) = 𝐹(1.4, −1) = ( ),
−1 0.2
and
[𝐽(𝑋0 )]−1 = [𝐽(1.4, −1)]−1 = (0.2885 −0.0962
)
0.1923 0.2692
Thus
𝑥1 1.4 0.2885 −0.0962 −0.04 1.4308
(𝑦 ) = ( ) − ( )( )=( )
1 −1 0.1923 0.2692 0.2 −1.0462
root to 2 d.p. is
𝑥1 = 1.43 and 𝑦1 = −1.05.
Further iterations can be done by updating 𝑋0 by 𝑋1.
𝑥2 + 𝑦 − 1
(e) Define the functions as 𝐹=[ 3 ]
𝑦 − 2𝑥 + 4
2𝑥 1
we have the Jacobian 𝐽=[ ].
−2 3𝑦 2
>> clear
>> F=@(x,y) [x.^2+y-1; y.^3-2*x+4];
>> J=@(x,y)[2*x, 1; -2, 3*y.^2];
>> x(1)=1.4;
2
Spring 2024-2025
>> y(1)=-1;
>> for n=1:50
Xn=[x(n); y(n)];
Fn=F(x(n), y(n));
Jn=J(x(n), y(n));
Xn1=Xn-Jn\Fn;
x(n+1)=Xn1(1);
y(n+1)=Xn1(2);
end
>> Roots=[x',y']
Roots =
1.4000 -1.0000
1.4308 -1.0462
1.4299 -1.0447
1.4299 -1.0447
1.4299 -1.0447
(f) >> clear
% Left_hand side of equations as vector
>> F=@(x) [x(1)^2+x(2)-1; x(2)^3-2*x(1)+4];
% Guess solution
>> x0=[1; -1];
% Solve the system
>> x=fsolve(F, x0)
x=
1.4299
-1.0447
4.2 Fixed Point Iteration Method for Nonlinear System
Consider the nonlinear system of two variables of the form
𝑓1 (𝑥, 𝑦) = 0 and𝑓2 (𝑥, 𝑦) = 0 (1)
By rearrangements the system can be expressed as
𝑥 = 𝑔1 (𝑥, 𝑦)and 𝑦 = 𝑔2 (𝑥, 𝑦) (2)
If there is a point (𝑝, 𝑞) such that
𝑝 = 𝑔1 (𝑝, 𝑞)and 𝑞 = 𝑔2 (𝑝, 𝑞)
then the point (𝑝, 𝑞) is a fixed point of the system and it is a root of the system.
With the arrangement (2), we may assume an iterative formula
𝑥𝑛+1 = 𝑔1(𝑥𝑛,𝑦𝑛) and 𝑦𝑛+1 = 𝑔2 (𝑥𝑛 , 𝑦𝑛 ) (3)
where𝑛 = 0, 1, 2, 3,∙ ∙ ∙ .
Without further details it may be mentioned that if (𝑥0 , 𝑦0 ) is close to the fixed point
(𝑝, 𝑞) and if
𝜕 𝜕
| 𝑔1 (𝑥0 , 𝑦0 )| + | 𝑔1 (𝑥0 , 𝑦0 )| < 1 (4𝑎)
𝜕𝑥 𝜕𝑦
and
𝜕 𝜕
| 𝑔2 (𝑥0 , 𝑦0 )| + | 𝑔2 (𝑥0 , 𝑦0 )| < 1 (4𝑏)
𝜕𝑥 𝜕𝑦
then the iterative formula (3) will converge to the fixed point (𝑝, 𝑞).
3
Spring 2024-2025
If the conditions (4) are not satisfied, the iterative process might diverge.
The above method can be extended for more than two variables.
Seidel Iteration
An improvement of the iteration process is the Seidel iteration similar to Gauss-Seidel
iteration for linear system. In Seidel iteration process the iterative formulas are
𝑥𝑛+1 = 𝑔1 (𝑥𝑛 , 𝑦𝑛 )
𝑦𝑛+1 = 𝑔2 (𝑥𝑛+1 , 𝑦𝑛 )
Example 4.2
Consider the system of equations 𝑦 = 𝑥 2 − 5𝑥 + 3
𝑥 2 + 4𝑦 2 = 4
(a) Plot the graphs of the above system using MATLAB function “ezplot(fun)”.
From your plot find the number of real roots and estimate the roots.
(b) A fixed point iteration formula is suggested to estimate root at (𝑥0 , 𝑦0 )
1
𝑥𝑛+1 = (𝑥𝑛2 − 𝑦𝑛 + 3)
5
1
𝑦𝑛+1 = (4 − 𝑥𝑛2 − 4𝑦𝑛2 + 10𝑦𝑛 )
10
(i) Verify whether the above iterative formula will converge to the root near
(𝑥0 , 𝑦0 ) = (0.4, 1.0). If converges, perform one iteration otherwise suggest
another fixed point iterative formula which converge to root.
(ii) Verify whether the above iterative formula will converge to the root near
(𝑥0 , 𝑦0 ) = (0.9, −1). If converges, perform one iteration otherwise suggest
another fixed point iterative formula which converge to the root.
(c) Write MATLAB commands to execute the iterative formula in b(ii) five times.
(a) MATLAB codes for ezplot
>> clear
>> figure
>> hold on
>> ezplot('y-x.^2+5*x-3')
>> ezplot('x.^2+4*y.^2-4')
>> hold off
From graph it can be seen that the system has two real roots near (0.4, 1.0)and
(0.9, −0.9).
(b) (i) From the iterative formula let us define
1
𝑔1 (𝑥, 𝑦) = 5 (𝑥 2 − 𝑦 + 3)
and
1
𝑔2 (𝑥, 𝑦) = 10 (4 − 𝑥 2 − 4𝑦 2 + 10𝑦)
Here
𝜕 2𝑥 𝜕 1
𝑔1 (𝑥, 𝑦) = , 𝑔1 (𝑥, 𝑦) = −
𝜕𝑥 5 𝜕𝑦 5
and
4
Spring 2024-2025
𝜕 2𝑥 𝜕 1
𝑔2 (𝑥, 𝑦) = − , 𝑔2 (𝑥, 𝑦) = (−8𝑦 + 10)
𝜕𝑥 10 𝜕𝑦 10
Near (0.4, 1. ), we have
𝜕𝑔1 𝜕𝑔1
| |+| | = |0.16| + |−0.2| = 0.36 < 1
𝜕𝑥 𝜕𝑦
and
𝜕𝑔2 𝜕𝑔2 −8(1) + 10
| |+| | = |−0.08| + | | = 0.28 < 1
𝜕𝑥 𝜕𝑦 10
The partial derivative tests are satisfied and hence the iterations converges to the root
near (0.4, 1.0).
We can use iterative formulas for the root near (0.4, 1.0) as
1
𝑥𝑛+1 = (𝑥𝑛2 − 𝑦𝑛 + 3)
5
1
𝑦𝑛+1 = (4 − 𝑥𝑛2 − 4𝑦𝑛2 + 10𝑦𝑛 )
10
In Seidel iteration the second iterative equation will be of the form
1 2
𝑦𝑛+1 = (4 − 𝑥𝑛+1 − 4𝑦𝑛2 + 10𝑦𝑛 )
10
Calculations with the above iterative formulas with 𝑥0 = 0.4 and 𝑦0 = 1.0 are given
below:
Jacobi iteration Seidel iteration
n xn yn n xn yn
0 0.4 1.0 0 0.4 1.0
1 0.432 0.984 1 0.432 0.9813
(ii) Near (0.9, -1.0) we have
𝜕𝑔1 𝜕𝑔1
| |+| | = |0.36| + |−0.2| = 0.56 < 1
𝜕𝑥 𝜕𝑦
𝑔1 (𝑥, 𝑦) satisfy the partial derivative test.
𝜕𝑔 𝜕𝑔 −8(−1)+10
But | 𝜕𝑥2 | + | 𝜕𝑦2 | = |−0.18| + | | = 1.98 > 1
10
Test fails and the convergence is not guaranteed.
Rearranging the second equation by
𝑥 2 + 4𝑦 2 − 4 + 10𝑦 = 10𝑦
we have
1 2
𝑦= (𝑥 + 4𝑦 2 − 4 + 10𝑦) .
10
Corresponding iterative function can be of the form
1 2
𝑔3 (𝑥, 𝑦) = (𝑥 + 4𝑦 2 − 4 + 10𝑦) .
10
and
𝜕 2𝑥 𝜕 1
𝑔3 (𝑥, 𝑦) = , 𝑔3 (𝑥, 𝑦) = (8𝑦 + 10)
𝜕𝑥 10 𝜕𝑦 10
At (0.9, −1),
𝜕𝑔3 𝜕𝑔3 1.8 8(−1) + 10
| |+| |=| |+| | = 0.38 < 1
𝜕𝑥 𝜕𝑦 10 10
5
Spring 2024-2025
𝑔1 𝑥, 𝑦) satisfies the partial derivative test.
Thus a fixed point iterative formula which converge to root near (0.9, −1) is
1
𝑥𝑛+1 = (𝑥𝑛2 − 𝑦𝑛 + 3)
5
1 2
𝑦𝑛+1 = (𝑥 + 4𝑦𝑛2 − 4 + 10𝑦𝑛 )
10 𝑛
(e) >> clear
>> x(1)=0.9;
>> y(1)=-1;
>> for n=1:5
x(n+1)=(x(n)^2-y(n)+3)/5;
y(n+1)=(x(n)^2+4*y(n)^2-4+10*y(n))/10;
end
>> Iterative_Roots =[x',y']
Iterative_Roots =
0.9000 -1.0000
0.9620 -0.9190
0.9689 -0.8886
0.9655 -0.8789
0.9622 -0.8767
0.9605 -0.8767
6
Spring 2024-2025
Exercise 4
1. Consider the following nonlinear systems:
(a) 𝑥 3 + 𝑦 2 = 2, 𝑥 2 − 𝑦 = 1 starting with (𝑥0 , 𝑦0 ) = (1.2, −0.5)
(b) 3𝑥 2 − 𝑦 2 = 0, 3𝑥𝑦 2 − 𝑥 3 − 1 = 0 starting with (𝑥0 , 𝑦0 ) = (1, 1).
(c) ln(𝑥 2 + 𝑦 2 ) + 𝑦 = 1, √𝑥 + 𝑥𝑦 2 = 0 starting with (𝑥0 , 𝑦0 ) = (2.4, −0.6).
1 1
(d) 4𝑥 2 − 20𝑥 + 4 𝑦 2 + 8 = 0, 2 𝑥𝑦 2 + 2𝑥 − 5𝑦 + 8 = 0 starting with (𝑥0 , 𝑦0 ) =
(0.5, 1.9).
(e) 𝑦 sin(𝑥𝑦) + 2 = 0, cos(𝑥𝑦) + 𝑥 − 𝑦 = 0 starting with (𝑥0 , 𝑦0 ) = (1, 2).
2
(f) tan−1 𝑥 + 𝑒 𝑦 = 5, cos 𝑥 2 + 𝑥𝑦 = 10 starting with (𝑥0 , 𝑦0 ) = (−1.6, −1.3).
Now using Newton-Raphson method answer the following questions. For each of the
above nonlinear system of equations:
(i) Find the Jacobian matrix 𝐽(𝑥, 𝑦) for the above system.
(ii)Evaluate the inverse of the Jacobian matrix at the given point.
(iii) Write down the iterative formula for Newton-Raphson method.
(iv) Estimate the root to 2 d.p. using the using the above iterative formula once starting
with the given point.
(v)Write down MATLAB commands to execute the iteration four times.
(vi) Use MATLAB function “fsolve (fun,x0)” to find the root.
.
2. Sketch the curves represented by the following non-linear equations,
2𝑒 𝑥 + 𝑦 = 0,
3𝑥 2 + 4𝑦 2 = 8.
a) Solve using Newton’s method with initial estimate (−1, −2).
b) Repeat (a) but with an initial estimate (−2,0).
c) Use graphical approach to validate the results in (a) and (b).
3. Sketch the curves represented by the following non-linear equations,
𝑥 2 − 𝑦 2 = 4,
𝑥 2 + 𝑦 2 = 16.
a) Solve using Newton’s method with initial estimate (3.0,2.5).
b) Repeat (a) but with an initial estimate (−3, −2.5).
c) Use graphical approach to validate the results in (a) and (b).
7
Spring 2024-2025
4. The system of non-linear equations
4𝑥 + ln(𝑥 2 + y) = 1 and 8𝑦 − ln(𝑥 + 𝑦 2 ) = 1
has a solution near (𝑥0 , 𝑦0 ) = (1, 2).
The following iterative formula is suggested to estimate the roots.
1
𝑥𝑛+1 = (1 − ln(𝑥𝑛2 + 𝑦𝑛 ))
4
1
𝑦𝑛+1 = 8 (1 + ln (𝑥𝑛 + 𝑦𝑛2 ))
i. Verify whether the above iterative formula will converge to the root near
(𝑥0 , 𝑦0 ) = (1, 2). If converges, perform ONE iteration otherwise write ‘Not
Convergent’.
ii. Write MATLAB commands to execute the iterative formula five times.
.
5. Determine the roots of the following simultaneous non linear equations using (a)
Newton-Raphson method, (b) fixed-point iteration method.
𝑥 = 𝑦 + 𝑥 2 − 0.5,
𝑦 = 𝑥 2 − 5𝑥𝑦
Employ initial guesses of 𝑥0 = 𝑦0 = 1.0 and discuss the results.
6. The following system has a root near (𝑥0 , 𝑦0 ) = (0.2, 0.3).
2𝑥 2 + 𝑦 2 − 15𝑥 + 2 = 0,
𝑥𝑦 2 + 𝑥 − 10𝑦 + 5 = 0.
Estimate the root correct to 3 decimal places using fixed point iterative method.
7. Problem # 3.13 Page- 90 [Ref. Numerical methods for Engineers and Scientists –
Amos Gilat, Vish Subramaniam].
8. Problem # 3.14 Page- 90 [Ref. Numerical methods for Engineers and Scientists –
Amos Gilat, Vish Subramaniam].
8
Spring 2024-2025
7. A planar, two-link robot arm is shown in the figure. The coordinate system 𝑥𝑦 is the
tool frame and attached to the end-effector. The coordinates of the end-effector relative to
the base frame are expressed as
𝑥 = 𝐿1 cos 𝜃1 + 𝐿2 cos(𝜃1 + 𝜃2 )
𝑦 = 𝐿1 sin 𝜃1 + 𝐿2 sin(𝜃1 + 𝜃2 )
Suppose the lengths, in consistent physical units, of the two links are 𝐿1 = 1 and 𝐿2 = 2,
and that 𝑥 = 2.5, 𝑦 = 1.4. Find the joint angles 𝜃1 and 𝜃2 (in radians) using Newton’s
method with initial estimate of (0.8,0.9)
8. The equation of the catenary curve and the ellipse, which are shown in the figure, are
given by:
1 𝑥 −𝑥
𝑓1 (𝑥, 𝑦) = 𝑦 − (𝑒 2 + 𝑒 2 ) = 0
2
𝑓2 (𝑥, 𝑦) = 9𝑥 2 + 252 − 225 = 0
Use Newton’s method to determine the point of intersection of
the curves that resides in the first quadrant of the coordinate
system. Assume the initial guess values as 𝑥0 = 2.5, 𝑦0 = 2.0.
[Ref. Numerical Methods for engineers and Scientists – Amos
Gilat, Vish Subramaniam, Example # 3.5, Page # 83]
9. Solve the following system of nonlinear equations:
−2𝑥 3 + 3𝑦 2 + 42 = 0
5𝑥 2 + 3𝑦 3 − 69 = 0
i. Use Newton’s Method. Start at 𝑥 = 1, 𝑦 = 1, (𝑥 = 3.5, 𝑦 = 2.5) and carry out
the two iterations.
ii. Use fixed-point iteration method. Use the iteration functions 𝑦 =
1/3 1/3
−5𝑥 2 +69 3𝑦 2 +42
( ) and 𝑥 = ( ) . Start 𝑥 = 1, 𝑦 = 1, and carry out the first six
3 2
iterations.
[Ref. Numerical Methods for engineers and Scientists – Amos Gilat, Vish
Subramaniam, Problem # 3.13, Page # 90]
10. Solve the following system of nonlinear equations:
𝑥 2 + 2𝑥 + 2𝑦 2 − 26 = 0
9
Spring 2024-2025
2𝑥 3 − 𝑦 2 + 4𝑦 − 19 = 0
i. Use Newton’s Method. Start at 𝑥 = 1, 𝑦 = 1, (1.5, 2) and carry out the two
iterations.
ii. Use fixed-point iteration method. Start 𝑥 = 1, 𝑦 = 1, (1.5, 2) and carry out
1 1
19+𝑦 2 −4𝑦 3 26−𝑥 2 −2𝑥 2
the first four iterations. ( 𝑥 = ( ) ,𝑦 = ( ))
2 2
[Ref. Numerical Methods for engineers and Scientists – Amos Gilat, Vish
Subramaniam, Problem # 3.14, Page # 90]
10
