NUMERICAL METHODS
WITH APPLICATIONS
(MEC500)
Tutorial 2
TAYLOR SERIES
EXPANSION
Question 1
Evaluate f(x) when x = 1.2 using TSE
f (x, y) 3x 15x 24
4 3
x0 = 1 base point
h = x – x0 step size
Solution f (x, y) 3x 15x 24 4 3
clear all
%to evaluate f(x) when x=1.2 using TSE
x = 1.2; % to evaluate f(1.2)
fe = -3*x^4 + 15*x^3 +24; % exact value
x0 = 1; %base point
h = x - x0; %step size
f0 = -3*x0^4 + 15*x0^3 +24; % f(x)
%derivatives of function, f(x)
fd(1) = -12*x0^3 + 45*x0^2; % f'(x)
fd(2) = -36*x0^2 + 90*x0; % f''(x)
fd(3) = -72*x0 + 90; % f'''(x)
fd(4) = -72; % f''''(x)
fd(5) = 0; % f'''''(x)
Solution (cont)
maxit = 5; % stopping criteria: max
% number of iterations
es = 0.05; % stopping criteria: desired
% tolerance, % es=0.05%.
fa(1,:) = f0; % zero order TSE
et(1,:) = abs((fe-fa(1,:))/fe)*100; % true error of zero order TSE
storedata(1,:) = [0 fa(1,:) et(1,:) inf]; % store result of zero order TSE
Solution (cont)
for i=1:maxit %stopping criteria: max no if
% iterations = 'maxit'
%TSE & truncating the remainder terms
fa(i+1,:) = fa(i) + fd(i)*(h^i)/factorial(i);
%calculating the true error, et
et(i+1,:) = abs((fe-fa(i+1,:))/fe)*100;
%calculating the approximate error, ea
ea(i+1,:) = abs((fa(i+1,:)-fa(i,:))/fa(i+1,:))*100;
%store all data in 'storedata'
storedata(i+1,:) = [i fa(i+1,:) et(i+1,:) ea(i+1,:)];
Solution (cont)
% stopping criteria: desired tolerance, es = 0.05%
if ea(i+1,:)<es
break; %stop iteration if desired tolerance is
% reached
end
end
storedata%#ok<NOPTS> %show data at the end of the iteration
ROOT FINDING
Bisection Method – example 1
• Find the root of equation below:
x
f (x) e x
using the bisection method with x1 = 0, x2 = 1, and ε =10-3.
xl xu
x new
2
x
f ( x) e x
Solution
Editor – new script (function)
function y = f(x)
y = exp(-x) - x;
Save as m.files f.m
Command window
>>clear
>> format long
eps_abs = 1e-5;
eps_step = 1e-5;
a = 0.0;
b = 1.0;
while (b - a >= eps_step || ( abs( f(a) ) >= eps_abs && abs( f(b) ) >= eps_abs ) )
c = (a + b)/2;
if ( f(c) == 0 )
break;
elseif ( f(a)*f(c) < 0 )
b = c;
else
a = c;
end
end
Command window - continue
>> [a b]
ans =
0.567138671875000 0.567146301269531
>> abs(f(a))
y=
7.237911846869061e-06
ans =
7.237911846869061e-06
>> abs(f(b))
y=
-4.718446080853589e-06
ans =
4.718446080853589e-06
Thus, the root is x = 0.56714
Bisection Method – example 2
• Find the root of equation below:
f ( x) x 4 x 10
3 2
using the bisection method with x1 = 1, x2 = 2, n=13
function[x e] = mybisect(f,a,b,n)
% function [x e] = mybisect(f,a,b,n)
Editor – new script
% Does n iterations of the bisection method for a function f
% Inputs: f --an inline function
% a,b --left and right edges of the interval
% n --the number of bisections to do.
% Outputs: x --the estimated solution of f(x) = 0
% e --an upperbound on the error
format long
c = f(a); d = f(b);
if c*d > 0.0
error('Function has same sign at both endpoints.')
end
disp(' x y')
for i = 1:n
x = (a + b)/2;
y = f(x);
disp([ x y])
if y == 0.0 % solved the equation exactly
e = 0;break % jumps out of the for loop
end
if c*y < 0
b=x;
else
a=x;
end
end Save as m.files mybisect.m
Command window:
To run in command window:
x y
1.500000000000000 2.375000000000000
>> f=@(x) x.^3+4*x.^2-10 1.250000000000000 -1.796875000000000
1.375000000000000 0.162109375000000
f= 1.312500000000000 -0.848388671875000
1.343750000000000 -0.350982666015625
@(x)x.^3+4*x.^2-10 1.359375000000000 -0.096408843994141
1.367187500000000 0.032355785369873
>> mybisect(f,1,2,13) 1.363281250000000 -0.032149970531464
1.365234375000000 0.000072024762630
1.364257812500000 -0.016046690754592
1.364746093750000 -0.007989262812771
1.364990234375000 -0.003959101522923
1.365112304687500 -0.001943659010067
ans =
1.365112304687500
False Position Method – example 1
Find the real root of the equation:
f ( x) x 3 2 x 2 5
by using false position method between x = 2 and 3.
Formula of the root:
f (xu )(xl xu )
xr xu
f (xl ) f (xu )
Let us find first approximation (x2) by taking Manually calculated
x0=2, x1=3,
X2 = x0 – {(x1 – x0)/[f(x1) – f(x0)]}*f(x0)
Therefore, x2 = 2 – [(3-2)/(16+1)]*(-1) = Therefore, x6=2.0939
2.0588 Here, f(x6)= -0.0075
Here, f(x2) = f(2.0588) = -0.3908 = negative Hence root lies between (2.0939,3)
Hence, roots lies between (2.0588,3)
Therefore, x7=2.0943
X3 = x1 – {(x2 – x1)/[f(x2) – f(x1)]}*f(x1) Here, f(x8)= -0.0027
Therefore, x3 = 2.0588 – [(3- Hence root lies between (2.0943,3)
2.0588)/(16+0.3908)]*(-0.3908) = 2.0813
Here, f(x3) = f(2.0813) = -0.1469 = negative Therefore, x9=2.0945
Hence roots lies between (2.0813,3). Here, f(x9)= -0.0010
Hence root lies between (2.0945,3)
Therefore, x4=2.0896
Here, f(x4)= -0.0547 Therefore, x10=2.0945
Hence root lies between (2.0896,3)
Now, x9 and x10 are equal, therefore we
Therefore, x5=2.0927 can stop here.
Here, f(x5)= -0.0202
Hence root lies between (2.0927,3) So our final answer is x=2.0945.
f=@(x) x^3-2*x-5; Editor – new script
a=2; b=3;
for i=1:10
x0=a; x1=b;
fprintf('\n Hence root lies between (%.4f,%.0f)',a,b)
x2(i)=x0-(x1-x0)/(f(x1)-f(x0))*f(x0);
if f(x2(i))>0
b=x2(i);
else a=x2(i);
end
fprintf('\n Therefore, x2=%.4f \n Here,
f(x2)=%.4f',x2(i),f(x2(i)))
p=x2(i);
end
for i=1:10
error(i)=p-x2(i);
end
Answer=p
plot (error)
grid on;
title('Plot of error');
xlabel('iterations');
ylabel('Error'); Save as m.files falseposition.m
Therefore, x2=2.0939
Command window Here, f(x2)=-0.0075
Hence root lies between (2.0939,3)
Therefore, x2=2.0943
>> falseposition Here, f(x2)=-0.0027
Hence root lies between (2.0943,3)
Hence root lies between (2.0000,3) Therefore, x2=2.0945
Therefore, x2=2.0588 Here, f(x2)=-0.0010
Here, f(x2)=-0.3908 Hence root lies between (2.0945,3)
Hence root lies between (2.0588,3) Therefore, x2=2.0945
Therefore, x2=2.0813 Here, f(x2)=-0.0004
Here, f(x2)=-0.1472 Hence root lies between (2.0945,3)
Hence root lies between (2.0813,3) Therefore, x2=2.0945
Therefore, x2=2.0896 Here, f(x2)=-0.0001
Here, f(x2)=-0.0547 Hence root lies between (2.0945,3)
Hence root lies between (2.0896,3) Therefore, x2=2.0945
Therefore, x2=2.0927 Here, f(x2)=-0.0001
Here, f(x2)=-0.0202
Hence root lies between (2.0927,3) Answer =
2.094546950874257
False Position & Fixed Point Iteration
Given
f ( x) x3 6.8x 2 13x 6
Find the roots using False Position and Fixed Point Iteration between
x = 0 and 5.
f (xu )(xl xu )
False Position: xr xu
f (xl ) f (xu )
Fixed Point Iteration:
1. Rearrange the function f(x) = 0 so that x is on the left side of the
equation.
2. Express x = g(x) in iterative formula to predict a new value of x.
Newton Raphson – example 1
• Find the root of equation
f ( x) cos x 3x 1
using the Newton-Raphson method with starting point x1 = 0.0
f(xi )
xi 1 = xi -
f '(xi )
f (xi 1 ) xi 1 xi
Editor – new script / command window
% Newton Raphson Method
clear all
close all
clc
% Change here for different functions
f=@(x) cos(x)-3*x+1
%this is the derivative of the above function
df=@(x) -sin(x)-3
% Change lower limit 'a' and upper limit 'b'
a=0; b=1;
x=a;
for i=1:1:100
x1=x-(f(x)/df(x));
x=x1;
end
Continue…
sol=x;
fprintf('Approximate Root is %.15f',sol)
a=0;b=1;
x=a;
er(5)=0;
for i=1:1:5
x1=x-(f(x)/df(x));
x=x1;
er(i)=x1-sol;
end
plot(er)
xlabel('Number of iterations')
ylabel('Error')
title('Error Vs. Number of iterations')
Save as m.files newtonraphson1.m
Command window
>> newtonraphson1
f=
@(x)cos(x)-3*x+1
df =
@(x)-sin(x)-3
Approximate Root is 0.607101648103123
Newton Raphson – example 2
• Find the root of equation
f ( x) x x 4 / 3
using the Newton-Raphson method with starting point x1 = 1.0 and
the convergence criterion ε =1x10-6, maximum number of
iterations=10
f(xi )
xi 1 = xi -
f '(xi )
f (xi 1 ) xi 1 xi
Editor – new script/command window
f = @(x) x + x^(4/3); % Equation definition
fp = @(x) 1 + 4/3*x^(1/3); % First-order derivative of f
x0 = 1; % Initial guess
N = 10; % Maximum number of iterations
tol = 1E-6; % Convergence tolerance
x = zeros(N + 1,1); % Preallocate solution vector where row
=> iteration
x(1) = x0; % Set initial guess
% Newton's Method algorithm
n = 2;
nfinal = N + 1; % Store final iteration if tol is reached
before N iterations
while (n <= N + 1)
fe = f(x(n - 1));
fpe = fp(x(n - 1));
x(n) = x(n - 1) - fe/fpe;
if (abs(fe) <= tol)
nfinal = n; % Store final iteration!break;
end
n = n + 1;
end
Continue…
% Plotevolution of the solution
figure('Color','White')
plot(0:nfinal - 1,x(1:nfinal),'o-')
title('Newton''s Method Solution: $f(x) = x +
x^{\frac{4}{3}}$','FontSize',20,'Interpreter','latex')
xlabel('Iteration','FontSize',16)
ylabel('$x$','FontSize',16,'Interpreter','latex')
Save as m.files newtonraphson2.m
Command window
>>newtonraphson2
>>sol=x
sol =
1.000000000000000
0.142857142857143
0.014668874758294
0.000902408208838
0.000025750234322
0.000000243864662
0.000000000503664
0.000000000000133
0.000000000000000
0.000000000000000
0.000000000000000
Secant method
• Find the root of equation f (x) e x x
using the secant method with x1 = 0, x2 = 1, and and the convergence
criterion ε =10-3.
f (x i )(xi x i 1 )
xi 1 x i
f (x i ) f (xi 1)
xi 1 xi f (xi 1) ,
Editor – new script
% Find the roots of using Secant Method
% func --> function is taken by user
% like @(x)(sin(x)) or @(x)(x^2 + x - 1) or any other function
% but use the same format i.e. use paranthesis as used above with '@' and 'x'
% maxerr --> Maximum Error in finding Root
maxiter = 100; % max number of iteration before get the answer
f = input('Enter Function in terms of x: ');
xn_2 = input('Ener Lower Limit: ');
xn_1 = input('Ener Upper Limit: ');
maxerr = input('Enter Maximum Error: ');
xn = (xn_2*f(xn_1) - xn_1*f(xn_2))/(f(xn_1) - f(xn_2));
disp('xn-2 f(xn-2) xn-1 f(xn-1) xn f(xn)');
disp(num2str([xn_2 f(xn_2) xn_1 f(xn_1) xn f(xn)],'%20.7f'));
flag = 1;
while abs(f(xn)) > maxerr
xn_2 = xn_1;
xn_1 = xn;
xn = (xn_2*f(xn_1) - xn_1*f(xn_2))/(f(xn_1) - f(xn_2));
disp(num2str([xn_2 f(xn_2) xn_1 f(xn_1) xn f(xn)],'%20.7f'));
flag = flag + 1;
if(flag == maxiter)
break;
end
end
if flag < maxiter
display(['Root is x = ' num2str(xn)]);
else
display('Root does not exist');
end
Save as m.files secant.m
Command window
>> clear
>> secant
Enter Function in terms of x:@(x)exp(-x)-x;
Enter Lower Limit:0
Enter Upper Limit:1
Enter Maximum Error:1E-6
xn-2 f(xn-2) xn-1 f(xn-1) xn f(xn)
0.0000000 1.0000000 1.0000000 -0.6321206 0.6126998 -0.0708139
1.0000000 -0.6321206 0.6126998 -0.0708139 0.5638384 0.0051824
0.6126998 -0.0708139 0.5638384 0.0051824 0.5671704 -0.0000424
0.5638384 0.0051824 0.5671704 -0.0000424 0.5671433 -0.0000000
Root is x = 0.56714
GAUSS METHOD
Gauss-Seidel Method
Find the roots of following simultaneous equations using the Gauss-Seidel
method:
20x + y – 2z = 17 …………………………...(1)
3x + 20y – z = – 18 …………………………...(1)
2x – 3y + 20z = 25 ……………………………(1)
Let us write the given equation in the form as:
x = (1/20)(17-y+2z)……………..(2)
y = (1/20)(-18-3x+z)…………….(2)
z = (1/20)(25-2x+3y)……………(2)
we start approximation by x0=y0=z0=0.
Substituting y=y0 and z=z0 in right hand side of first of equation (2)
X1 = 17/20=0.85
In second of equation second put x=x1 and z=z0
Y1 = (1/20)(-18-3*0.85) = -1.0275
Put x=x1, y=y1 in third of equation (2)
Z1 = (1/20)(25-2*0.85+3*-1.0275) = 1.011
Similarly we get,
x2=1.002; y2 =-0.9998; z2 =0.9998
x3 =1.0000; y3 =-1.0000; z3 =1.0000
x4 =1.0000; y4 =-1.0000; z4 =1.0000
If you observe above two sets of roots, they are almost same.
Hence the roots of given simultaneous equations using Gauss-Seidel Method are:
x=1; y = -1; z=1
% Find the roots of given
simultaneous equations using Gauss
Seidel Methods: >> x
f1=@(x,y,z) (1/20)*(17-y+2*z); x =
f2=@(x,y,z) (1/20)*(-18-3*x*z);
f3=@(x,y,z) (1/20)*(25-2*x+3*y); 1.001702192138035
x=0; y=0; z=0; >> y
for i=1:5 y =
f1(x,y,z);
f2(x,y,z); -1.049122726589640
f3(x,y,z);
x=f1(x,y,z); >> z
y=f2(x,y,z);
z=f3(x,y,z); z =
end
0.992461371797750
Gauss Method
Given:
4 1 1 0 x1 275
1 4 0 1x 225
2
1 0 4 1x3 75
0 1 1 4 x 25
4
Find the x1, x2, x3 and x4 using Gauss Elimination and Gauss Seidel
Method.
