0% found this document useful (0 votes)
117 views

Numerical Analysis Assignment

The document contains code for implementing numerical methods like Gauss-Seidel method, Secant method, and Newton-Raphson method to find the roots of nonlinear equations. It includes code snippets submitted by multiple students to solve equations using these three methods. The code takes the coefficient matrix and source vector as input, calculates the roots iteratively, and outputs the solution vector.
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
0% found this document useful (0 votes)
117 views

Numerical Analysis Assignment

The document contains code for implementing numerical methods like Gauss-Seidel method, Secant method, and Newton-Raphson method to find the roots of nonlinear equations. It includes code snippets submitted by multiple students to solve equations using these three methods. The code takes the coefficient matrix and source vector as input, calculates the roots iteratively, and outputs the solution vector.
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
You are on page 1/ 14

Numerical Analysis Lab

Submitted to: Ma’am Kanwal


Submitted by: Muhammad Shahood Jamal (2019-CH-15)
M Jawad (2019-CH-17)
Faraz Afzal (2019-CH-55)
Naseer Khan (2018-CH-77)
Fahad Ali (2016-CH-76)

Department of Chemical Engineering


Page 1 of 14
Submitted By: M Jawad (2019-CH-17)
Submitted To: Mam Kanwal

Newton Method:
clear all
clc
N = 100;
y =@(x) x^3+2*x^2+10; %fd = sym('x^3+2*x^2+10’); f2 =diff(fd)
yd =@(x) 3*x^2+4*x;
x = 0;
X(1) = x;
err =0.01;
for i = 2:N
x = x - y(x)/yd(x);
X(i) = x; I = i-1 %it will update the new value of x
Err = abs(X(i) - X(i-1)); if Err < err, break; end
disp(['The Root is: ' num2str(x) ' ,with accuracy: '
num2str(Err) ' ,No iterations: ' num2str(I)])
end

Secant Method:
clear all
clc
%s= sin(x)-4*x+6 first point of interval= 0 ,second point of
interval=-0.01 Error= 0.001
a=input('Enter function:','s');
f=inline(a)

x(1)=input('Enter first point of guess interval: ');


x(2)=input('Enter second point of guess interval: ');
n=input('Enter allowed Error in calculation: ');
iteration=0;

for i=3:1000
x(i) = x(i-1) - (f(x(i-1)))*((x(i-1) - x(i-2))/(f(x(i-1)) -
f(x(i-2))));
iteration=iteration+1;
if abs((x(i)-x(i-1))/x(i))*100<n
root=x(i)
iteration=iteration
break

Page 2 of 14
end
end

Seidel Method:
clear all
clc
% 7*x1-0.2*x2-0.3*x3=21
% 0.3*x1+10*x2-0.5*x3=31.5
% 0.4*x1-0.5*x2+12*x3=51.5
%Initial values: x2=0, x3=0
%iterate until the value of error for x1 is less than 0.01%
%solution:
%x1=(21+0.2*x2+0.3*x3)/7
%x2=(31.5-0.3*x1+0.5*x3)/10
%x3=(51.5-0.4*x1+0.5*x2)/12

clear all
clc; format('long','g')
i=1;
x2(i)=0; x3(i)=0;
error_x1(1)=9999;
while error_x1(i) >=0.01;
x1(i+1)=(21+0.2*x2(i)+0.3*x3(i))/7;
x2(i+1)=(31.5-0.3*x1(i+1)+0.5*x3(i))/10;
x3(i+1)=(51.5-0.4*x1(i+1)+0.5*x2(i+1))/12;
error_x1(i+1)=abs((x1(i+1)-x1(i))/x1(i+1))*100; %absolute
error=abs(approximate value-exact value)/approximate value*100
error_x2(i+1)=abs((x2(i+1)-x2(i))/x2(i+1))*100;
error_x3(i+1)=abs((x3(i+1)-x3(i))/x3(i+1))*100;
i=i+1;
end
disp(' x1 error(%)');
disp([x1' ,error_x1']);
disp(' x2 error(%)');

disp([x2' ,error_x2'])
disp(' x3 error(%)');

disp([x3' ,error_x3'])

Page 3 of 14
Submitted By: M Shahood Jamal (2019-CH-15)
Submitted To: Mam Kanwal

GAUSS SEIDAL METHOD

n=input('Enter number of equations, n: ');


A = zeros(n,n+1);
x1 = zeros(n);
tol = input('Enter the tolerance, tol: ');
m = input('Enter maximum number of iterations, m: ');
A=[4 2 3 8; 3 -5 2 -14; -2 3 8 27];
x1=[0 0 0];

k = 1;
while k <= m
err = 0;
for i = 1 : n
s = 0;
for j = 1 : n
s = s-A(i,j)*x1(j);
end
s = (s+A(i,n+1))/A(i,i);
if abs(s) > err
err = abs(s);
end
x1(i) = x1(i) + s;
end

if err <= tol


break;
else
k = k+1;
end
end
fprintf('Solution vector after %d iterations is :\n', k-1);
for i = 1 : n
fprintf(' %11.8f \n', x1(i));
end

Page 4 of 14
SECANT METHOD
d=input('Enter function:','s');
f=inline(d)
x(1)=input('Enter first point of guess interval: ');
x(2)=input('Enter second point of guess interval: ');
i=input('Enter allowed Error in calculation: ');
iteration=0;
for n=3:100
x(n)=(x(n-2)(f(x(n-1)))-x(n-1)(f(x(n-2))))/(f(x(n-1))-f(x(n-2)));
iteration=iteration+1;
iteration=iteration
x(n)
if abs((x(n)-x(n-1))/x(n-1))<i
break
end
end

NEWTON RAPHSONS METHOD


x = 0.2;
fx = 3*x+sin(x)-exp(1)^x;
f_prime = 3+ cos(x)-x*exp(1)^x;
x1 = 2;
if ((fx ~= 0) && (f_prime~= 0))
while abs(x-x1)>= 0.000001
x1 = x - fx/f_prime;
x = x1;
end

end
fprintf('The root of the equation is : %f\n',x1);

Page 5 of 14
Submitted By: Faraz Afzal (2019-CH-55)
Submitted To: Mam Kanwal

GAUSS SEIDAL METHOD

A=input('Enter a coefficient matrix A:');


B=input('Enter source vector B:');
P=input('Enter initial guess vector:');
n=input('Enter number of iterations');
e=input('Enter your tolerance');
N=length(B);
X=zeros(N,1);
Y=X+1; %for stopping criteria
j=1
while abs(Y-X)>e
Y=X
for i=1:N
X(i)=(B(i)/A(i,i))- (A(i,[1:i-1,i+1:N])*P([1:i-1,i+1:N]))/A(i,i);
P(i)=X(i);
end
fprintf('Iteration no %d\n',j)
X
j=j+1
end

SECANT METHOD
f= @(x) x^3-6*(x^2)+11*x-6.1;
x0= 4;
x1= 7;

Page 6 of 14
i=0.001;
iteration=0;
for n=1:100
x2 = (x0*f(x1) - x1*f(x0)) / (f(x1) - f(x0));
iteration=iteration+1;
iteration=iteration
x2
if abs((x2-x1)/x1)<i
break
end
x0=x1;
x1=x2;
end
end

NEWTON RAPHSONS METHOD


fx =@(x) 3*x+sin(x)-exp(1).^x;
f_prime = diffs(fx,x);

while abs(x-x1)>= 0.000001


x1 = x - fx(x)/f_prime;
x = x1;
f_prime = diffs(fx,x);
end

function dv = diffs(fx,x)
h = 0.05;
for i = 10
r = (fx(x+h)+fx(x))/h;

Page 7 of 14
h = h/2;
end
dv = r;
end

Submitted By: Naseer (2018-CH-77)


Submitted To: Mam Kanwal

GAUSS SEIDAL METHOD


disp('Give the input to solve the set of equations AX=B')
a=input('Input the square matrix A : \n');
b=input('Input the column matrix B : \n');
m=length(a);
%z is a two dimensional array in which row corresponds to values of X in a
%specific iteration and the column corresponds to values of specific
%element of X in different iterations
c=0;
%random assignment
e=1;
%'e' represents the maximum error d=0;
%random assignment
for u=1:m
x(u)=b(u,1)/a(u,u);
z(1,u)=0;
%initializing the values for matrix X(x1;x2;...xm) end l=2;
%'l' represents the iteration no.
%loop for finding the convergence factor (C.F) for r = 1:m for s = 1:m if r~=s
p(r)=abs(a(r,s)/a(r,r))+d;%p(r) is the C.F for equation no. r d=p(r); end end d=0; end if
min(p)>=1 %atleast one equation must satisfy the condition p<1 fprintf('Roots will not
converge for this set of equations') else while(e>=1e-4) j1=1;%while calculating elements

Page 8 of 14
in first column we consider only the old values of X for i1=2:m
q(j1)=(a(j1,i1)/a(j1,j1))*z(l-1,i1)+c;
c=q(j1);
end
c=0;
z(l,j1)=x(j1)-q(j1);
%elements of z in the iteration no. l x(j1)=z(l,j1);
for u=1:m
x(u)=b(u,1)/a(u,u);
z(1,u)=0;
end
for j1=2:m-1
%for intermediate columns between 1 and m, we use the updated values of X
for i1=1:j1-1
q(j1)=(a(j1,i1)/a(j1,j1))*z(l,i1)+c;
c=q(j1);
end
for i1=j1+1:m q(j1)=(a(j1,i1)/a(j1,j1))*z(l-1,i1)+c; c=q(j1);
end
c=0; z(l,j1)=x(j1)-q(j1); x(j1)=z(l,j1);
for u=1:m x(u)=b(u,1)/a(u,u); z(1,u)=0;
end
j1=m;
%for the last column, we use only the updated values of X
for i1=1:m-1
q(j1)=(a(j1,i1)/a(j1,j1))*z(l,i1)+c; c=q(j1);
end
c=0;
z(l,j1)=x(j1)-q(j1);
Page 9 of 14
for v=1:m
t=abs(z(l,v)-z(l-1,v));
%calculates the error
end
e=max(t);
%evaluates the maximum error out of errors of all elements of X l=l+1;
%iteration no. gets updated
for i=1:m
X(1,i)=z(l-1,i);
%the final solution X
end
end
%loop to show iteration number along with the values of z
for i=1:l-1
for j=1:m
w(i,j+1)=z(i,j);
end
w(i,1)=i;
end
disp(' It. no. x1 x2 x3 x4 ')
disp(w)
disp('The final solution is ')
disp(X)
fprintf('The total number of iterations is%d',l-1)

SECANT METHOD
d=input('Enter function:','s');
f=inline(d)
x(1)=input('Enter first point of guess interval: ');

Page 10 of 14
x(2)=input('Enter second point of guess interval: ');
i=input('Enter allowed Error in calculation: ');
iteration=0;
error=9999;
n=3;
while error>i;
x(n) = (x(n-2)(f(x(n-1)))-x(n-1)(f(x(n-2))))/(f(x(n-1))-f(x(n-2)));
iteration=iteration+1;
iteration=iteration
x(n)
error= abs((x(n)-x(n-1))/x(n-1));
n=n+1;
end

NEWTON RAPHSONS METHOD


% Clearing Screen
clc
% Setting x as symbolic variable
syms x;

% Input Section
y = input('Enter non-linear equations: ');
a = input('Enter initial guess: ');
e = input('Tolerable error: ');
N = input('Enter maximum number of steps: ');
% Initializing step counter
step = 1;
g = diff(y,x);
fa = eval(subs(y,x,a));

Page 11 of 14
while abs(fa)> e
fa = eval(subs(y,x,a));
ga = eval(subs(g,x,a));
if ga == 0
disp('Division by zero.');
break;
end

b = a - fa/ga;
fprintf('step=%d\ta=%f\tf(a)=%f\n',step,a,fa);
a = b;

if step>N
disp('Not convergent');
break;
end
step = step + 1;
end

fprintf('Root is %f\n', a);

Submitted By: Fahad Ali (2016-CH-76)


Submitted To: Mam Kanwal

GAUSS SEIDAL METHOD


A=input('Enter a coefficient matrix A:');
B=input('Enter source vector B:');
P=input('Enter initial guess vector:');
n=input('Enter number of iterations');
e=input('Enter your tolerance');

Page 12 of 14
N=length(B);
X=zeros(N,1);
Y=zeros(N,1); %for stopping criteria
for j=1:n
for i=1:N
X(i)=(B(i)/A(i,i))- (A(i,[1:i-1,i+1:N])*P([1:i-1,i+1:N]))/A(i,i);
P(i)=X(i);
end
fprintf('Iteration no %d\n',j)
X
if abs(Y-X)<e
break
end
Y=X
end

SECANT METHOD
d=input('Enter function:','s');
f=inline(d)
x(1)=input('Enter first point of guess interval: ');
x(2)=input('Enter second point of guess interval: ');
i=input('Enter allowed Error in calculation: ');
iteration=0;
error=9999;
for n=3:100
if error > i
x(n)=(x(n-2)(f(x(n-1))) - x(n-1)(f(x(n-2)))) / (f(x(n-1)) - f(x(n-2)));
iteration=iteration+1;
iteration=iteration

Page 13 of 14
x(n)
error = abs((x(n)-x(n-1))/x(n-1));
else
break
end
end

NEWTON RAPHSONS METHOD


fx =@(x) 3*x+sin(x)-exp(1).^x;
f_prime = diffs(fx,x);

while abs(x-x1)>= 0.000001


x1 = x - fx(x)/f_prime;
x = x1;
f_prime = diffs(fx,x);
end

function dv = diffs(fx,x)
h = 0.05;
for i = 10
r = (fx(x+h)+fx(x))/h;
h = h/2;
end
dv = r;
end

Page 14 of 14

You might also like