MOHAMED SATHAK ENGINEERING COLLEGE

KILAKARAI - 623 806

DEPARTMENT OF CHEMICAL ENGINEERING

CX5111 COMPUTATIONAL PROGRAMME IN CHEMICAL ENGINEERING

LABORATORY

M.TECH - CHEMICAL ENGINEERING

2017 REGULATION
 Contents

S.No DATE EXPERIMENT SIGNATURE

1. BISECTION METHOD

2. FALSE POSITION METHOD

3. NEWTON RAPHSON METHOD

4. SIMPSON’S 1/3RD RULE

5. SIMPSON’S 3/8TH RULE

6. JACOBBI’S ITERATION METHOD

7. HEAT EXCHANGER OPTIMATION

8. TRAPEZOIDAL RULE

9. GAUSS ELIMINATION METHOD

10. GAUSS JORDAN METHOD

11. LU FACTORIZATION

12. SIMULATION OF CSTR

13. ADAM BASHFORTH METHOD

14. EULER’S SERIES

 Bisection method
Ex.no: Date:

Aim:

To find out the roots by using bisection method

Algorithm:

1. Input function and limits.

2. Repeat steps 3 and 4 100
times. 3. x=(a+b)/2
4. If f(x0)<0, a=x else b=x
5. Display x
6. Repeat steps 7 and 8 10 times.
7. Error = x-(a+b)/2
8. Store error values in array
9. Plot error
10. STOP.

Example

Find the roots of x^3-x-1 using bisection method and plot the error

MATLAB Code:

f=@(x) x^3-x-
1; a=1;
b=2;
for
i=1:100
c=(a+b)/2
; if f(c)>0
b=c;
else
a=c;
end
end
a=1; b=2;
p=c; for
i=1:100
c=(a+b)/2;
er(i)=f(c)-
f(p); if f(c)>0
b=c;
else
a=c;
end
end
fprintf('Root of given equation is %f',c) plot(er);
title('Plot of error')
 xlabel('Number of iterations')
ylabel('Error')
grid on;:

Result:

MATLAB programme code successfully verified in Bisection method

 False Position Method
Ex.no: Date:

Aim:

To find out the real roots of equations of f(X)=0

False Position Method MATLAB Code:

Find the roots of x^3-2*x-5=0 by using false position method.

clc;

f=@(x) x^3-2*x-5; 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');
Result:

MATLAB programme code Successfully verified

 Newton Raphson Method
Ex.no: Date:

Aim:

To find out the real roots of equations of f(X)=0

MATLAB Program:

% 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

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')

Result:

MATLAB programme code Successfully verified

 Simpson’s 3/8th Rule MATLAB Program
Ex.no: Date:

AIM:

To findout the roots from the equation of f(x)=0

MATLAB Code:

Example:

%Question: Evaluate the integral x^4 within limits

-3 to 3

clc;

clear all;

close all;

f=@(x)x^4; %Change here for different

function a=-3;b=3; %Given limits

n=b-a; %Number of

intervals h=(b-a)/n;

p=0;

for i=a:b

p=p+1; x(p)=i;

y(p)=i^4; %Change here for different

function end

l=length(x)

;x

y answer=(3*h/8)*((y(1)+y(l))+3*(y(2)+y(3)+y(5)+y(6))

+2*(y(4)))
Result:

MATLAB programme code Successfully verified

 JACOBBI ITERATION METHOD

Ex.no: Date:

AIM:

To find out the roots of the equation using Jacobbi iteration method in Matlab

MATLAB Code:

%Jacobi’s Iterative Method

% find the roots of given simultaneous equations using Jacobi’s

method: clc;

f1=0(x,y,z) (1/20)*(17-y+2*z);

f2=0(x,y,z) (1/20)*(-16-3*x+z);

f3=0(x,y,z) (1/20)*(25-2*x+3*z);

x=0; y=0; z=0;

for i=1:5

f1(x,y,z);

f2(x,y,z);

f3(x,y,z);

a=f1 (x,y,z) ;

b=f2(x,y,z);

c=f3(x,y,z);

x=a; y=b;

z=c; end

z
Result:

MATLAB programme code Successfully verified

 HEAT EXCHANGER-OPTIMATION

Ex.no: Date:

Aim:

To Design optimization of Heat Exchanger

MATLAB Code:

function [f1 x] = ftrialr2 (x)

% Ds = x(1);
% do = x(2);
% B = x(3);
a1 = 8000;
a2 = 259.2;
a3 = 0.93;
n = 1; %number of passes
eff = 0.7; %Pump
effficiency Nt = 60; %no.
of tubes
% tube side parameters
mt = 18.80;
Tci = 37.8;
Tco = 76.7;
rhot = 995.0;
mut = 0.00358;
%viscosity muwt =
0.00213;
cpt = 2.05;
kt = 0.13;
Rft = 0.00061;
st = 1.25*x(2);
vt = (mt/(rhot*0.8*x(2).^2*pi/4))*n/Nt; %velocity of fluid on
tube side di = 0.8*x(2);
Ret = rhot*vt*0.8*x(2)/mut; %Reynold's number
tube side ft = 0.079/(Ret^0.25); %Darcy friction
factor
Prt = mut*cpt/kt; %Prandtl number tube side
% shell side
parameters ms =
5.52;
Thi = 199.0;
Tho = 93.3;
rhos = 850.0;
cps = 2.47;
mews = 0.0004;
mewws = 0.00036;
ks = 0.13;
Rfs = 0.00061;
ce = 0.12; %energy cost
H = 7000; %annual operating time in hours
de = 4*(0.43*st.^2 - (0.5*pi*x(2).^2/4))/(0.5*pi*x(2));
%equivalent dia As = x(1)*x(3)*(1-x(2))/(1.25*x(2)); %shell
side crooss section area
vs = ms/(rhos*As); %velocity of fluid on shell side
Res = ms*de/(As*mews); %shell side Reynold's
number Prs = mews*cps/ks; %shell side Prandtl's
number
%Thermal Calculations
hs = 0.36*(kt/de)*(Res^0.55)*(Prs^(1/3))*(mews/mewws)^0.14;
%shell side heat transfer coefficient
R = (Thi - Tho)/(Tco - Tci); %correction
coefficient P = (Tco - Tci)/(Thi - Tci);
%efficiency

F = sqrt((R.^2+1)/(R.^2-1)); %correction factor

LMTD = ((Thi - Tco) - (Tho -
Tci))/(log((Thi-Tco)/(Tho-Tci))); Q = ms*cps*(Thi -
Tho);
L = 4.5;
syms ht U A L
[ht, U, A, L] = solve((ht == kt/di*(3.657 +
(0.0677*(Ret*Prt*((di/L).^1.33)).^1/3))/(1+0.1*Prt*(Ret*(di/L)).^0.3)),
U == 1/((1/hs)+Rfs+(x(2)/di)*(Rft+(1/ht))), A == Q/(U*F*LMTD),
L == A/(pi*x(2)*Nt));
if Ret < 2300
ht = kt/di*(3.657 +
(0.0677*(Ret*Prt*((di/L).^1.33).^(1/3)))/(1+0.1*Prt*(Ret*(di/L)).^
0.3)); else if Ret > 2300 && Ret < 10000
ht = kt/di*((1+di/L).^0.67*(ft/8)*(Ret -
1000)*Prt/(1+12.7*((ft/8).^0.5)*((Prt.^(2/3))-1)));
else if Ret > 10000
ht =
0.027*kt/do*(Ret.^0.8)*(Prt.^(1/3)).*(mewt/mewwt).^0.14;
end
en
d
en
d
U = 1/((1/hs)+Rfs+
(x(2)/di)*(Rft+(1/ht))); A =
Q/(U*F*LMTD);
L = A/(pi*x(2)*Nt);
%Pressure drop
%tube side pressure drop
p = 4; %p = 2.5; %constant
pt = 0.5*rhot*vt.^2*(L*ft/(0.8*x(2))+p)*n;
%shell side pressure
drop b0 = 0.72;
fs = 2*b0*Res;
ps = fs*(rhos*vs.^2/2)*(L/x(2))*(x(1)/de);
f1 = (a1 + a2*A^a3)+ (ce*H*(mt*pt/rhot +
ms*ps/rhos)/eff); end
function [hist,searchdir] =
runfmincon global hist
hist.x = [];
hist.fval = [];
searchdir = [];
x0 = [0.3 0.02 0.2]
lb = [0.2 0.015 0.05 3];
ub = [2 0.051 0.5 6];
options = optimset('fmincon');
options = optimset(options,'Display','Iter-detailed');
options = optimset('PlotFcns',@optimplotfvalcust)
figure
options =
optimset('PlotFcns',@optimplotxcust) options
= optimset(options,'outputfcn',@outfun);
nonlcon = [];
 [x,f1val] = fmincon(@ftrialr2,x0,[],[],[],
[],lb,ub,nonlcon,options)

end

Result:

MATLAB programme code Successfully verified

 SIMPSON’S 1/3RD RULE

Ex.no: Date:

AIM:

To Findout the roots of given equation of f(x)=0

MATLAB Code:

%Question: Evaluate the integral X^4 within limits 3 to

-3 clc;

clear all;

close all;

f=@(x)x^4; %Change here for different

function a=-3;b=3; %Given limits

n=b-a; %Number of

intervals h=(b-a)/n;

p=0;
for

i=a:b

p=p+1

x(p)=i

y(p)=i^4; %Change here for different

function end

l=length(x);

x
y

answer=(h/3)*((y(1)+y(l))+2*(y(3)+y(5))+4*(y(2)+y(4)+y(6)))
Result:

MATLAB programme code Successfully verified

 TRAPEZOIDAL RULE

Ex.no: Date:

AIM:

To Evaluate the roots of the given equations

MATLAB Code:

%Question: Evaluate the integral X^4 within limits 3 to

-3 clc;

clear all;

close all;

f=@(x)x^4; %Change here for different

function a=-3;b=3; %Given limits

n=b-a; %Number of

intervals h=(b-a)/n;

p=0;
for i=a:b

p=p+1;

x(p)=i;

y(p)=i^4; %Change here for different

function end

l=length(x);

answer=(h/2)*((y(1)+y(l))+2*(sum(y)-y(1)-y(l)))
Result:

MATLAB programme code Successfully verified

 GAUSS ELIMINATION METHOD

Ex.no: Date:

AIM:

To findout the roots of the given equation

f(x)=0 MATLAB Code:

[m,n]=size(a);

for j=1:m-1
for z=2:m
if a(j,j)==0
t=a(j,:);a(j,:)=a(z,:); a(z,:)=t;
end
end
for i=j+1:m
a(i,:)=a(i,:)-a(j,:)*(a(i,j)/a(j,j));
end
end
x=zeros(1,m);

for s=m:-1:1 c=0;

for k=2:m
c=c+a(s,k)*x(k);
end
x(s)=(a(s,n)-c)/a(s,s);
end
disp('Gauss elimination method:'); a
x'

Result:

MATLAB programme code Successfully verified

 GAUSS –JORDAN METHOD

Ex.no: Date:

AIM:
To findout the roots of the given equation

MATLAB Code:

% Gauss-Jordan method
[m,n]=size(a);

for j=1:m-1
for z=2:m
if a(j,j)==0
t=a(1,:);a(1,:)=a(z,:);
a(z,:)=t;
End
End
for i=j+1:m
a(i,:)=a(i,:)-a(j,:)*(a(i,j)/a(j,j));
End
end
for j=m:-1:2
for i=j-1:-1:1
a(i,:)=a(i,:)-a(j,:)*(a(i,j)/a(j,j));
end
end
for s=1:m
a(s,:)=a(s,:)/a(s,s);

x(s)=a(s,n);

end

disp('Gauss-Jordan method:');

a
x'
Result:

MATLAB programme code Successfully verified

 LU FACTORIZATION
Ex.no: Date:

AIM:

To find out the roots of given equation by triangular matrix

MATLABCode:

function [L,U] =

my_lu(A) n = size(A,

1);

I = eye(n); L = I;

U = A; for

k=1:n-1

L(k+1:n,k) = U(k+1:n,k)/U(k,k);

%multipliers for j = k+1:n

U(j,k:n) = U(j,k:n)-L(j,k)*U(k,k:n);

%rows end

end

function [L,U] = my_lu2(A)

n = size(A, 1); I = eye(n); L = I; U

= A; for k=1:n-1

L(k+1:n,k) = U(k+1:n,k)/U(k,k); %multipliers

U(k+1:n,k:n) = U(k+1:n,k:n)-L(k+1:n,k)*U(k,k:n); %block

function [L,A]=LU_factor(A,n)
% LU factorization of an n by n matrix A
% using Gauss elimination without pivoting
% LU_factor.m
% A is factored as A = L*U
% Output:
% L is lower triangular with the main diagonal part = 1s.
% U is upper triangular and is stored in the original mtx A
% and must be zeroed out to get U
% K. Ming Leung, 01/26/03

L=eye(n);
for k=1:n
if (A(k,k) == 0) Error('Pivoting is
needed!'); end
L(k+1:n,k)=A(k+1:n,k)/A(k,k);
for j=k+1:n
A(j,:)=A(j,:)-L(j,k)*A(k,:);
end
end
Result:

MATLAB programme code Successfully verified

 SIMULATION OF CSTR

Ex.no: Date:

AIM:

To modelling reactor CSTR by simulink

MATLAB Code :

function f = cstrrxns(C)

load cstrrxns.mat

% The data stored in cstrrxns.mat are:

%k1 = 1.0; k2 = 0.2; k3 = 0.05; k4 = 0.4; Cao = 1; t = 2;

%F1 = -Ca + Cao [-k1Ca - k2Ca^3/2 + k3Cc^2]t = 0

%F2 = -Cb + [2k1Ca - k4Cb^2]t = 0
%F3 = -Cc + [k2Ca^3/2 - k3Cc^2 + k4Cb^2]t = 0
%F4 = -Cd + [k4Cb^2]t = 0

f(1) = -C(1) + Cao + (-k1*C(1) - k2* C(1)*sqrt(C(1)) + k3*C(3)*C(3))*t; f(2) =

-C(2) + (2*k1*C(1) - k4*C(2)*C(2))*t;
f(3) = -C(3) + (k2*(C(1)^(3/2)) - k3*C(3)*C(3) + k4*C(2)*C(2))*t; f(4) =
-C(4)+ ( k4*C(2)*C(2) )*t;

% Initial values for iteration:

%C0 = [1 0 0 0]'

\
Result:

MATLAB programme code Successfully verified

 ADAM BASHFORTH METHOD
Ex.no: Date:

AIM:

To run the predictor corrector - Adam Bashforth method using C Program.

C-Program:

C PROGRAM FOR ADAM BASHFORTH METHOD

#include<stdio.h>
#include<conio.h>
#include<math.h>
#define MAX 20

float equation(float x,float y)

{
return(1+(y*y));

}
void main()
{
FILE *fp;
int i=0,count=-1;
float lower,upper,h,y1,xvalue[MAX],yvalue[MAX],result;
float function[MAX],search,final,temp;
fp=fopen("admbsh.dat","w");
clrscr();
printf("ADAM-BASHFORTH METHOD ");
fprintf(fp,"ADAM-BASHFORTH METHOD ");
printf("\n");
fprintf(fp,"\n");
printf("\nEnter the lower bound of x : ");
fprintf(fp,"\nEnter the lower bound of x : ");
scanf("%f",&lower);
fprintf(fp,"%f",lower);
printf("\nEnter the upper bound of x : ");
fprintf(fp,"\nEnter the upper bound of x : ");
scanf("%f",&upper);
fprintf(fp,"%f",upper);
printf("\nEnter the value of y(lower) : ");
fprintf(fp,"\nEnter the value of y(lower) : ");
scanf("%f",&y1);
fprintf(fp,"%f",y1);
printf("\nEnter the value of h : ");
fprintf(fp,"\nEnter the value of h : ");
scanf("%f",&h);
fprintf(fp,"%f",h);
printf("\nEnter the value of x for which you want to find y :");
fprintf(fp,"\nEnter the value of x for which you want to find y :");
scanf("%f",&search);
fprintf(fp,"%f",search);
xvalue[i]=lower;
 yvalue[i]=y1;
for(i=0;xvalue[i]<=upper;i++)
{
xvalue[i+1]=xvalue[i]+h;
}
for(i=0;xvalue[i]<=upper;i++)
{
result=equation(xvalue[i],yvalue[i]);
yvalue[i+1]=yvalue[i]+(h*result);
printf("\n\n");
fprintf(fp,"\n\n");
printf("The table is ");
fprintf(fp,"The table is ");
printf("\n\n");
fprintf(fp,"\n\n");
printf(" i x y f(x,y) ");
fprintf(fp," i x y f(x,y) ");
printf("\n\n");
fprintf(fp,"\n\n");
for(i=0;xvalue[i]<=upper;i++)
{
function[i]=equation(xvalue[i],yvalue[i]);
printf(" %d. %.4f %.4f %.4f ",i,xvalue[i],yvalue[i],function[i]);
fprintf(fp," %d. %.4f %.4f %.4f ",i,xvalue[i],yvalue[i],function[i]);
count=count+1;
printf("\n");
fprintf(fp,"\n");
}
yvalue[search]=yvalue[count]+(h/24)*((-9*function[count-3])+(37*function[count-2])-
(59*function[count-1])+(55*function[count]));
final=equation(search,yvalue[search]);
yvalue[search]=yvalue[count]+(h/24)*((function[count-2])-(5*function[count-1])
+(19*function[count])+(9*final));
printf("\n\n");
fprintf(fp,"\n\n");
printf("Approximate value is : %.4f ",yvalue[search]);
fprintf(fp,"Approximate value is : %.4f ",yvalue[search]);
fclose(fp);
getch();
}
Result:

C- Program for Adam Bashforth predictor corrector is successfully verified

 ADAM BASHFORTH METHOD
Ex.no: Date:

AIM:

To run the Euler’s series using C Program.

C-Program:

C-PROGRAM FOR EULER’S SERIES

#include<stdio.h>
float fun(float x,float y)
{
    float f;
    f=x+y;
    return f;
}
main()
{
    float a,b,x,y,h,t,k;
    printf("\nEnter x0,y0,h,xn: ");
    scanf("%f%f%f%f",&a,&b,&h,&t);
    x=a;
    y=b;
    printf("\n  x\t  y\n");
    while(x<=t)
    {
        k=h*fun(x,y);
        y=y+k;
        x=x+h;
        printf("%0.3f\t%0.3f\n",x,y);
    }
}

Result:

C- Program for Euler’s series is successfully verified.