Newton’s forward interpolation
function fx= forward(x,y,b)
x=input('x in metrix form');
y=input('y in metrix form');
b=input('point to be interpolated');
n=length(x);
%coad for table construction
for j=1:n-1
if j==1
for i=1:n-1
del(i,1)=y(i+1)-y(i);
end
else
for i=1:n-j
del(i,j)=del(i+1,j-1)-del(i,j-1);
end
end
end
del
h=x(2)-x(1);
p=(b-x(1))/h;
%cod for newton's forward interpolation formula
for j=1:n-1
if j==1
c(1)=p*del(1,j);
else
t=p;
for k=1:j-1
t=t*(p-k);
end
c(j)=(t*del(1,j))/factorial(j);
end
end
sum(c);
ans=y(1)+sum(c)
end
Backword
function fx= backward(x,y,b)
x=input('x in metrix form');
y=input('y in metrix form');
b=input('point to be interpolated');
n=length(x);
%coad for table construction
for j=1:n-1
if j==1
for i=1:n-1
del(i,1)=y(i+1)-y(i);
end
else
for i=1:n-j
del(i,j)=del(i+1,j-1)-del(i,j-1);
end
end
end
del
h=x(2)-x(1);
p=(b-x(n))/h;
% coding for newton's backword interpolation formula
for j=1:n-1
if j==1
c(1)=p*del(n-1,j);
else
t=p;
for k=1:j-1
t=t*(p+k);
end
c(j)=(t*del(n-j,j))/factorial(j)
end
end
sum(c);
ans=y(n)+sum(c)
Newton’s divided difference interpolation
function fp=devided_interpolation(x,y,p)
clc;
x=input('enter the value of x in vector form ');
y=input('enter the value of y in vector form ');
p=input('enter the point to be interpolated ');
n=length(x);
for i=1:n-1
d(i,1)=(y(i+1)-y(i))/(x(i+1)-x(i));
end
for j=2:n-1
for i=1:n-j
d(i,j)=(d(i+1,j-1)-d(i,j-1))/(x(i+j)-x(i));
end
end
d
%coefficient loop
z(1)=1;
c(1)=y(1);
for j=2:n
z(j)=(p-x(j-1))*z(j-1);
a(j)=d(1,j-1);
c(j)=a(j)*z(j);
end
fp=sum(c);
end
Weddles rule
a=input('enter lower limit ');
b=input('enter upper limit ');
n=input('no of subdivision');
h=(b-a)/(n-1);
i=0;
for p=1:n
x=a+(p-1)*h;
y=1/(1+(x*x));
if(p==1)||(p==n)
co=(3/10);
else if (mod(p,6)==4)
co=(9/5);
else if (mod(p,6)==0)||(mod(p,6)==2)
co=(3/2);
else co=(3/10);
end
end
end
i=i+co*y*h;
end
i
simpson’s one third
clc;
a=input('the value of lower limit=');
b=input('the value of upper limit=');
n=input('no of subdivision');
h=(b-a)/n;
int=0;
for p=0:n
x = a + p*h;
y = 1/(1 + x*x);
if (p==0)||(p==n)
coeff=(1/3);
else if (mod(p,2)==0)
coeff=(2/3);
else
coeff=(4/3);
end
end
int=int+coeff*y*h;
end
int
simpson 3/8 rule
clc;
a=input('the value of lower limit=');
b=input('the value of upper limit=');
n=input('no of subdivision');
h=(b-a)/n;
I=0;
for p=0:n
x = a + p*h;
y = 1/(1 + x*x);
if (p==0)||(p==n)
coeff=(3/8);
else if (mod(p,3)==0)
coeff=(3/4);
else
coeff=(9/8);
end
end
I=I+coeff*y*h;
end
I
Trapezoidal rule
clc;
clear all;
a=input('the value of lower limit=');
b=input('the value of upper limit=');
n=input('number of values=');
h=(b-a)/(n-1);
I=0;
for p=1:n
x(p) = a + (p-1)*h;
y(p) = 1/(1 + x(p)*x(p));
if (p==1)||(p==n)
coeff=0.5;
else
coeff=1;
end
I=I+coeff*y(p)*h;
plot(x,y,'*')
end
I
Simpson’s 1/3 with data entry
clc;
a=input('the value of lower limit=');
b=input('the value of upper limit=');
table=xlsread('123.xlsx');
x=table(:,1);
y=table(:,2);
%x=input('enter the value of x in vector form');
% y=input('enter the value of y in vector form');
n=length(x);
h=(b-a)/(n-1);
int=0;
for i=1:n
if (i==1)||(i==n)
coeff(i)=(1/3);
else if (mod(i,2)==0)
coeff(i)=(4/3);
else
coeff(i)=(2/3);
end
end
int=int+coeff(i)*y(i)*h
end
int
Runge-kutta method
clc;
clear all;
x0= input ('inital value of x = ');
y0= input ('initial value of y = ');
xn= input ('final value of x = ');
h= input ('increment in x = ');
n= (xn-x0)/h;
f=@(x,y)(y-2*x/y);
for i=1:n
x=x0+(i-1)*h;
k1=h*f(x,y0);
k2=h*f(x+h/2,y0+k1/2);
k3=h*f(x+h/2,y0+k2/2); for loop
k4=h*f(x+h,y0+k3);
y1=y0 + (1/6)*(k1+2*k2+2*k3+k4);
y0=y1;
end
y1
clc;
clear all;
x0=input('initial value of x ');
y0=input('initial value of y ');
xf=input('final value of x');
n=input('no of steps ');
h=(xf-x0)/n;
ff=@(x,y)(y-2*x/y);
while(x0<xf)
k1=h*feval(ff,x0,y0);
k2=h*feval(ff,x0+h/2,y0+k1/2);
k3=h*feval(ff,x0+h/2,y0+k2/2);
k4=h*feval(ff,x0+h,y0+k3); while loop
y1=y0+1/6*(k1+2*(k2+k3)+k4);
y0=y1;
x1=x0+h;
x0=x1;
end
y1
Modified euler’s
clc;
clear all;
x0=input('initial value of x ');
y0=input('initial value of y ');
xf=input('final value of x');
n=input('no of steps ');
h=(xf-x0)/n;
ff=@(x,y) (x+y);
while(x0<xf)
yp=y0+h*feval(ff,x0,y0);
x1=x0+h;
yc=y0+(h/2*(feval(ff,x0,y0)+feval(ff,x1,yp)));
while(abs(yc-yp)>.0000005)
yp=yc;
yc=y0+(h/2*(feval(ff,x0,y0)+feval(ff,x1,yp)));
end
y0=yc;
x0=x1;
end
yc
Euler’s method
clc;
clear all;
x0=input('initial value of x ');
y0=input('initial value of y ');
xf=input('final value of x');
h=input ('step size = ');
% ff=input('enter function= ');
n=(xf-x0)/h;
ff=@(x,y) (x+y);
while(x0<=xf)
y1=y0+h*feval(ff,x0,y0);
x1=x0+h
y0=y1;
x0=x1
end
y1
Adam’s method
clc;
clear all;
x0=input('input initial value of x');
y0=input('input initial value of y');
xn=input('input final value of x');
n = input('number of step');
h=(xn-x0)/n;
z(1)=y0;
p(1)=x0;
fx = @(x,y) (y-(2*x)/y);
for i=2:4
k1 = h*feval(fx,x0,y0);
k2 = h*feval(fx,x0+h/2,y0+k1/2);
k3 = h*feval(fx,x0+h/2,y0+k2/2);
k4 = h*feval(fx,x0+h,y0+k3);
y1= y0+(k1+2*k2+2*k3+k4)/6;
z(i)=y1;
x0 = x0+h;
p(i)=x0;
y0=y1;
end
for i=5:11
p(i)=p(i-1)+h;
zp=z(i-1)+(h/24)*(55*fx(p(i-1),z(i-1))-59*fx(p(i-2),z(i-
2))+37*fx(p(i-3),z(i-3))...
-9*fx(p(i-4),z(i-4)));
z(i)=z(i-1)+(h/24)*(9*fx(p(i),zp)+19*fx(p(i-1),z(i-1))-
5*fx(p(i-2),z(i-2))...
+fx(p(i-3),z(i-3)));
while(abs(z(i)-zp)<0.000005)
zp=z(i);
z(i)=z(i-1)+(h/24)*(9*fx(p(i),z(i))+19*fx(p(i-1),z(i-
1))-5*fx(p(i-2),z(i-2))...
+fx(p(i-3),z(i-3)));
end
end
z(i)
Bisection (for loop)
function [x, e] = bisection(f,a,b,n)
f=@(x) x.^3+4*x.^2-10
a=input('enter the first end point: ');
b=input('enter the second end point: ');
n=input('enter the number of iterations: ');
format long
c=f(a);
d=f(b);
if(c*d)>0
error('function has same sign for both the end points')
end
disp(' x y ')
for i=1:n
x=(a+b)/2;
y=f(x);
disp([x y])
if y==0
e=0;
break
end
if(c*y)<0
b=x;
else
a=x;
end
end
plot(x,y)
Bisection (while loop)
function [x,e] = bisection(f,a,b,n)
z = input('Enter the function in x variable: ','s');
a = input('Enter the first end point: ');
b = input('Enter the second end point: ') ;
f = inline(z); %to use the expression as a function for a variable
t = input('Enter the tolerance to terminate iterations ');
format long
i=0;
c =f(a);
d =f(b);
if (c*d)>0
error('Function has same sign for both the end points')
end
disp (' x y ') % to differentiate
bw the values of x and y
while abs(b-a)>=t
i= i+1;
x = (a+b)/2;
y = f(x);
disp([x y]) %to show the iterations
if y == 0
e = 0;
break
end
if c*y < 0
b=x;
else
a=x;
end
end
str = ['The required root of the equation is: ', num2str(x), '']
NR method(for loop)
function x = mynewton (f,f1 ,x0 ,n,tol)
f =@(x) x.*tan(x)-1; % Inputs : f -- the function
f1 = @(x) x.*sec(x).*sec(x)+tan(x); % f1 -- it 's derivative
x0 = input('enter the value of x0:') % x0 -- starting guess , a number
tol = input ('enter the tolarance tol:') % tol -- desired tolerance , n =
input('enter the value of n:') % n - no.of iteration
% Output : x -- the approximate solution
x = x0; % set x equal to the initial guess x0
disp (' x y')
for i = 1:n % Do n times
x = x - f(x)/ f1(x);% Newton 's formula
y=f(x);
disp([ x y ])
plot(x,y,'*')
hold on
end
r = abs(f(x));
if r > tol
warning ('The desired accuracy was not attained ')
end
NR Method(while loop)
function [x] = mynewton(f,f1,x0,n,tol)
f=@(x)x^3-5;
f1=@(x)3*x^2;
x0=input('Enter the value of x0:')
tol=input('Enter the value of tolerance tol:')
n=input('Enter the value of n:')
format long
x=1:100
x=x0
y=f(x)
disp(' x y ')
while abs(y)>tol
x=x-f(x)/f1(x);
y=f(x);
disp([x y])
end
end
Secant method (for Loop)
function x = mysecant (f,x0 ,x1 ,n)
f=@(x)x.^3+4*x.^2-10
x0 = input('enter the value of x0:')
x1 = input ('enter the value of x1:')
n = input('enter the value of n:')
% Solves f(x) = 0 by doing n steps of the secant method
% starting with x0 and x1.
% Inputs : f -- the function
% x0 -- starting guess , a number
% x1 -- second starting guess
% n -- the number of steps to do
% Output : x -- the approximate solution
y0 = f(x0 );
y1 = f(x1 );
disp (' x y')
for i = 1:n % Do n times
x = x1 - (x1 -x0 )* y1 /(y1 -y0); % secant formula .
y=f(x);% y value at the new approximate solution .
disp([ x y ])
% Move numbers to get ready for the next step
x0=x1;
y0=y1;
x1=x;
y1=y;
end
end
Secant Method (while loop)
function[x]=mysecant(f,a,b,n,tol)
f=@(x)x*tan(x)-1
a=input('Enter a:')
b=input('Enter b:')
n=input('Enetr n:')
tol=input('Enter tol:')
x=a
y=f(x)
disp(' x y ')
while abs(y)>tol
x=a-((f(a)/(f(b)-f(a))*(b-a)));
y=f(x);
disp([x y])
end
end
PDE (liebman method)
clear all;
clc;
xdim=100;
ydim=100;
F_now=zeros(xdim+1,ydim+1);
F_prev=zeros(xdim+1,ydim+1);
i=1:1:xdim+1;
F_now(i,ydim+1)=45*((i-1)/xdim).*(1-((i-1)/xdim));
iter=0;
error=max(max(abs(F_now-F_prev)));
while(error>0.001)
iter=iter+1;
for i=2:1:xdim
for j=2:1:ydim
F_now(i,j)=F_now(i-1,j)+F_now(i+1,j)+F_now(i,j-1)+F_now(i,j+1)/4;
end
end
error=max(max(abs(F_now-F_prev)));
F_prev=F_now;
imagesc(F_now);colorbar;
title (['Fluid fLow Distribution in term of
velocity',int2str(xdim),'x',int2str(ydim),'Grid of Petroleum Reservoir at
iteration no',int2str(iter)],'Color','k');
getframe;
end
figure;
[ex,ey]=gradient(F_now);
quiver(-ex,-ey);
