Numerical Analysis Labs
Lab 1.1 Area by Double Integration
f = @(x, y) x.^2 + y.^2
a_x = 0;
b_x = 2;
a_y = @(x) x ;
b_y = @(x) 2 * x;
x_limits = [a_x, b_x];
y_limits = @(x) [a_y(x), b_y(x)];
area = integral2(f, a_x, b_x, a_y, b_y);
disp([‘The area under the curve is: ‘, num2str(area)]);
LAB 1.2 Volume
f = @(x, y, z) 1-(x.^2)/4-(y.^2)/9-(z.^2)/16
a_x = 0;
b_x = 2;
a_y = @(x) 0*x;
b_y = @(x) 3*sqrt(1-x.^2/4);
a_z = @(x,y) 0*x+0*y;
b_z = @(x,y) 4*sqrt(1-x.^2/4-y.^2/9);
x_limits = [a_x,b_x];
y_limits = @(x) [ a_y(x),b_y(x)];
z_limits = @(x,y) [ a_z(x,y),b_z(x,y),];
Volume = 8*integral3(f, a_x, b_x, a_y, b_y, a_z, b_z);
disp([‘The Volume under the curve is: ‘, num2str(Volume)]);
Lab 2 Improper Integrals
f = @(x) 1./(1+x.^2);
a = -Inf;
b = Inf;
integral_value = integral(f, a, b);
fprintf(‘The value of the improper integral is: %.6f\n’, integral_value);
Lab 3.1 Gradient
syms x y z
F=input(‘Enter a scalar field to obatain the gradient:’)
disp (‘Result:’)
gF=[diff(F,x),diff(F,y),diff(F,z)]
gradient=subs(gF,[x,y,z],[1,1,1])
Lab 3.2 Divergent and curl
syms x y z t
F=[(z^2+2*x*y) (x^2+2*y*z) (y^2+2*z*x)]
a(x,y,z)=F(1)
b(x,y,z)=F(2)
c(x,y,z)=F(3)
disp(‘The vector is:’);
F=[a,b,c];
disp([‘F’,’ = ‘])
disp(F)
dx=diff(a,x);
dy=diff(b,y);
dz=diff(c,z);
disp(‘The divergence is:’);
div=dx+dy+dz
colmx=diff(c,y)-diff(b,z);
colmy=diff(a,z)-diff(c,x);
colmz=diff(b,x)-diff(a,y);
disp(‘The curl is:’);
curl=[colmx,colmy,colmz]
curlF=subs(curl,[x,y,z],[1,-1,1])
Lab 6.1 Regula Falsi Method
syms x
f = input(‘Enter your function:’ )
a =input(‘Enter left side of your root guess: ‘ )
b = input(‘Enter right side of your root guess: ‘ )
n = input(‘Enter the number of iteration you want:’)
e = input (‘Enter your desired tolerance:’)
if f(a)*f(b)
for i = 1:n
c = (a*f(b)- b*f(a))/(f(b)- f(a));
fprintf(‘c%d=%.4f\n’,i,c)
%if abs(f(c))
%break
%end
if f(a)*f(c)
b=c;
elseif f(b)*f(c)
a=c;
end
end
fprintf(‘\nRoot is: %.4f\n’, c)
else
disp(‘No root between given bracket’)
end
Lab 6.2 Newton Raphson Method
syms x
f(x)=input(‘enter your function:’)
df=diff(f)
e=input(‘enter tolerance:’)
x0=input(‘enter initial value:’)
n=input(‘enter no of iterations:’)
if df(x0)~=0
for i=1:n
x1=x0-f(x0)/df(x0);
fprintf(‘x%d=%.4f\n’,i,x1)
if abs(x1-x0)
break
end
if df(x1)==0
disp(‘N-R method failed’)
end
x0=x1;
end
fprintf(‘\n root is:%.4f\n’,x1)
else
disp(‘N-R method failed’)
end
lab 7 Newtons forward and backward formula
x = [100 150 200 250 300 350 400]
y = [10.63 13.03 15.04 16.81 18.42 19.90 21.27]
n = length(x)
d = zeros(n-1);
h = x(2)-x(1)
xf = 218
p=(xf-x(1))/h
for k = 2:n % Calculation of first forward differences
d(k-1,1)=y(k)-y(k-1);%(y2-y1)
end
for r = 2:n-1% Calculation of second and rest forward differences
for k = 1:n-r
d(k,r) = d(k+1, r-1)-d(k, r-1);
end
end
disp(‘The difference table is :’)
d
s=y(1); t=p;
for r=1:n-1%Calculation
s=s+t*d(1,r);
t=(p-r)/(r+1)*t;
end
fprintf(‘The required value is f(%.1f)=%.4f’, xf,s);
Lab 8.1 Trapezoidal
syms x
f(x)=1/(1+x)
a=input(‘Enter lower limit a: ‘);
b=input(‘Enter upper limit b: ‘);
n=input(‘Enter the number of sub-intervals n: ‘);
h=(b-a)/n;
i = 1:1:n-1
xi=f(a+(i*h))
I = (h/2)*(f(a)+f(b)+2*sum(xi));
fprintf(‘The approximation of above integral is:%f’,I);
Lab 8.2 Simsons 1/3rd rule
syms x
f(x)=1/(1+x)
a=input(‘Enter lower limit a: ‘);
b=input(‘Enter upper limit b: ‘);
n=input(‘Enter the number of sub-intervals n: ‘);
if rem(n,2)==1
fprintf(‘\n Enter valid n!!!’)
n=input(‘\n Enter n as multiple of 2: ‘)
end
h=(b-a)/n;
i = 1:1:n-1 X=f(a+i*h)
even=sum(X(2:2:n-1))
odd=sum(X(1:2:n-1))
I = (h/3)*((f(a)+f(b)+4*odd +2*even))
fprintf(‘The approximation of above integral is:%f ‘,I);
Lab 8.3 Simpsons 3/8th rule
syms x
a = 0 % Lower Limit
b = 6 % Upper Limit
n = 6 % Number of Segments
f(x) = 1/(1+x.^2) % Declare the function
h = (b – a)/n % h is the segment size
xi=a:h:b
yi=f(xi)
i = 1:1:n-1;
S=f(a+i*h);
I=3:3:n-1;
S3=sum(S(I)); % compute sum
S(I)=[]; % Delete the values at I position
SO=sum(S); % compute sum
I = (3*h/8)*((f(a)+f(b)+3*SO + 2*S3));
fprintf(‘The approximation of above integral is %f\n’,I);
Lab 9.1 Taylors series
syms x y(x)
y1 = x-y^2 %define the function
x0 = 0
y0 = 1
X = 0.2
h = 0.1
y2 = diff(y1)
y3 = diff(y2)
y4 = diff(y3)
f1=diff(y,x,1)
f2=diff(y,x,2)
f3=diff(y,x,3)
y10 = subs(y1,{x,y(x)},[x0,y0])
y20 = subs(y2,{x,y(x),f1},[x0,y0,y10])
y30 = subs(y3,{x,y(x),f1,f2},[x0,y0,y10,y20])
y40 = subs(y4,{x,y(x),f1,f2,f3},[x0,y0,y10,y20,y30])
for i=x0+h:h:X
y = y0 + (x-x0)*y10 + (((x-x0)^2)/2)*y20 + (((x-x0)^3)/6)*y30 + (((x-x0)^4)/24)*y40
Y = subs(y,x,i);
fprintf(‘Value of y at x=%0.1f is %.4f\n’,i,Y)
end
Lab 9.2 Modified Eulers Method
f=@(x,y)x+y
x0=0
y0=1
xn=0.4
h=0.1
while x0
yE=y0+h*f(x0,y0);
x1=x0+h;
yM=y0+h/2*(f(x0,y0)+f(x1,yE));
err=abs(yM-yE);
x0=x1;
y0=yM;
table(x0,yM,y0)
end
fprintf(‘\n The value of y at x=%0.1f is y=%.4f’,xn,yM);
Lab 10.1 RK Method
f=@(x,y)y^2-x
x=input(‘Enter the initial value of x:’);
y=input(‘Enter the initial value of y(x):’);
h=input(‘Enter the step size(h):’);
X=input(‘Enter X at which Y is required:’);
for i=x+h:h:X
K1=double(h.*f(x,y))
K2=double(h.*f(x+h/2,y+K1/2))
K3=double(h.*f(x+h/2,y+K2/2))
K4=double(h.*f(x+h,y+K3))
K=(K1+2.*K2+2.*K3+K4)./6;
x=x+h;
y=y+K;
fprintf(‘Value of y at x=%0.2f is %f\n’,x,y)
end
Lab 10.2 Milnes Predictor method
syms x y
f(x,y) = 1+x*y^2
x0=0
x1=0.1
x2=0.2
x3=0.3
x4=0.4
h=0.1
y0=1
y1=1.105
y2=1.223
y3=1.355
y11=double(f(x1,y1))
y21=double(f(x2,y2))
y31=double(f(x3,y3))
y4=y0+4*h/3*(2*y11-y21+2*y31)
y41=double(f(x4,y4))
M=y2+(h/3)*(y21+4*y31+y41)
fprintf(‘The value is =%.4f\n’,M)