1. Iterative Calculation
a. Program for iterative calculation.
a. Program for iterative calculation.
clear;
n=3;
es=0.5*(10^(2-n));
x=0.5;
f(1)=1; //first estimate f=e^x=1
ft=1.648721; //true value of e^0.5=f
et(1)=(ft-f(1))*100/ft;
ea(1)=100;
i=2;
while ea(i-1)>=es
f(i)=f(i-1)+(x^(i-1)) //(factorial(i-1));
et(i)=(ft-f(i))*100/ft;
ea(i)=(f(i)-f(i-1))*100/f(i);
i=i+1;
end
for j=1:i-1
disp(ea(j),"approximate estimate of error(%)="et(j),"true % frelative error=",f(j),"result=",j,term number=)
disp(".....................................")
end
b. Program to calculate the roots of a quadratic equation using the formula.
responce=1;
while responce==1
a=input("input the value of a:")
b=input("input the value of b:")
c=input("input the value of c:")
if a==0 then
if b~=0 then
r1=-c/b;
disp(r1,"the root:")
else disp("trival solution")
end
else
discr=b^2-4*a*c;
if discr>=0 then
r1=(-b+sqrt(discr))/(2*a);
r2=(-b-sqrt(discr))/(2*a);
disp(r2,"and",r1,"the roots are:")
else
r1=-b/(2*a)
r2=r1;
i1=sqrt(abs(discr))/(2*a);
i2=i1;
disp(r2+i2*sqrt(-1),r1+i1*sqrt(-1),"the roots are:")
end
end
responce=input("do you want to continue (press 1 for yes and 2 if responce=2 then exit;)")
end
end
end
c. Program to evaluate 𝑒𝑥 using infinite series.
function y=f(x)
y=exp(x)
endfunction
sum=1;
test=0;
i=0;
term=1;
x=input("input value of x:")
while sum~=test
disp(sum,"sum:",term,"term:",i,"i:")
disp("...................")
i=i+1;
term=term*x/i;
test=sum;
sum=sum+term;
end
disp(f(x),"exact value")
2. Solution of algebraic and transcendental equations:
a. Program to solve algebraic and transcendental equation by bisection method.
deff('y=f(x)','y=x^3-x-1');
x1=1,x2=2; // f(1) is negative and f(2) is positive
d=0.0001; // for accuracy of root
c=1;
printf('successive approximations \tx1 \t \tx2 \t \tm \t \f f(m) \n')
while abs (x1-x2)>d
m=(x1+x2)/2;
printf('\t%f \t%f \t%f \t%f \n',x1,x2,m,f(m));
if f(m)*f(x1)>0
x1=m;
else
x2=m;
end
c=c+1; // to count number of iterations
end
printf(' the solution of equation after %i iteration is %g',c,m)
b. Program to solve algebraic and transcendental equation by false position method.
deff('y=f(x)','y=x^3-2*x-5');
a=2,b=3; //f(2) is negative and f(3) is positive
d=0.00001; // for accuracy of root
printf('successive iterations \ta \tb \t f(a) \t f(b) \tx1 \n')
for i=1:25
x1=b*f(a)/(f(a)-f(b))+a*f(b)/(f(b)-f(a))
if (f(a)*f(x1))>0
b=x1;
else
a=x1;
end
if abs(f(x1))<d
break
end
printf('\t%f %f %f %f %f\n',a,b,f(a),f(b),x1);
end
printf(' the root of equation %f',x1);
c. Program to solve algebraic and transcendental equation by Secant method.
deff('y=f(x)','y=x^3-2*x-5');
x1=2,x2=3
n=1;
c=0;
printf('successive iterations\tx1\tx2\tx3\tf(x3)\n')
while n==1
x3=(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1));
printf('\t%f\t%f\t%f\t%f\n',x1,x2,x3,f(x3));
if f(x3)*f(x1)>0 then
x2=x3;
else
x1=x3;
end
if abs(f(x3))<0.000001 then
break;
end
c=c+1;
end
printf('the root of the equation after %i iteration is:%f',c,x3)
d. Program to solve algebraic and transcendental equation by Newton Raphson method.
deff('y=f(x)','y=sin(x)-x/2');
deff('y1=f1(x)','y1=cos(x)-1/2');
x0=2,
d=0.0001;
c=0;n=1;
printf('successive iterations \tx0 \tf(x0) \tf1(x0) \n');
while n==1
x2=x0;
x1=x0-(f(x0)/f1(x0));
x0=x1;
printf('\t%f \t%f \t%f \n',x2,f(x1),f(x1));
c=c+1;
if abs (f(x0))<d
break;
end
end
printf(' the root of %i iteration is : %f',c,x0);
3. Interpolation
a. Program for Newton’s forward interpolation.
clc;clear;close
x=[1 3 5 7];
y=[24 120 336 720];
h=2
c=1;
for i=1:3
d1(c)=y(i+1)-y(i);
c=c+1;
end
c=1;
for i=1:2
d2(c)=d1(i+1)-d1(i);
c=c+1;
end
c=1;
for i=1:1
d3(c)=d2(i+1)-d2(i);
c=c+1;
end
d=[d1(1) d2(1) d3(1)];
x0=8;
pp=1;
y_x=y(1);
p=(x0-1)/2;
for i=1:3
pp=1;
for j=1:i
pp=pp*(p-(j-1))
end
y_x=y_x+(pp*d(i))/factorial(i);
end
printf('value of function at %f is:%f',lx0,y_x);
b. Program for Newton’s backward interpolation.
clc;clear;close
x=[15 20 25 30 35 40];
y=[0.2588190 0.3420201 0.4226183 0.5 0.5735764 0.6427876];
h=5
c=1;
for i=1:5
d1(c)=y(i+1)-y(i);
c=c+1;
end
c=1;
for i=1:4
d2(c)=d1(i+1)-d1(i);
c=c+1
end
c=1;
for i=1:3
d3(c)=d2(i+1)-d2(i);
c=c+1;
end
c=1;
for i=1:2
d4(c)=d3(i+1)-d3(i);
c=c+1;
end
c=1;
for i=1:1
d5(c)=d4(i+1)-d4(i);
c=c+1 ;
end
d=[d1(5) d2(4) d3(3) d2(2) d1(1)];
x0=38;
pp=1;
y_x=y(6);
p=(x0-x(6))/h;
for i=1:5
pp=1;
for j=1:i
pp=pp*(p-(j-1))
end
y_x=y_x+(pp*d(i))/factorial(i);
end
printf('value of function at %f is:%f',x0,y_x);
c. Program for Lagrange’s interpolation.
y=[4 12 19];
x=[1 3 4];
y_x=7;
Y_X=0;
poly(0,'y');
for i=1:3
p=x(i);
for j=1:3
if i~=j then
p=p*((y_x-y(j))/(y(i)-y(j)))
end
end
Y_X=Y_X+p;
end
disp(Y_X,'Y_X=');
4. Solving linear system of equations by iterative methods
a. Program for solving linear system of equations using Gauss Jordan method.
A=[2,1,1,10;3,2,3,18;1,4,9,16];
for i=1:3
j=i
while (A(i,i)==0&j<=3)
for k=1:4
B(1,k)=A(j+1,k)
A(j+1,k)=A(i,k)
A(i,k)=B(1,k)
end
disp(A);
j=j+1;
end
disp(A);
for k=4:-1:i
A(i,k)=A(i,k)/A(i,i)
end
disp(A)
for k=1:3
if (k~=i) then
l=A(k,i)/A(i,i)
for m=i:4
A(k,m)=A(k,m)-l*A(i,m)
end
end
end
disp(A)
end
for i=1:3
printf('\nx(%i)=%g\n',i,A(i,4))
end
b. Program for solving linear system of equations using Gauss Seidel method.
function [X]=gaussseidel(A, n, N, X, b)
L=A;
U=A;
D=A;
for i=1:1:n
for j=1:1:n
if j>i then L(i,j)=0;
D(i,j)=0;
end
if i>j then U(i,j)=0
D(i,j)=0;
end
if i==j then Ls(i,j)=0;
end
end
end
for k=1:1:N
X=(D+L)^-1*(-U*X+b);
disp(X)
end
endfunction
A=[2 -1 0; -1 2 -1; 0 -1 2]
b=[7;1;1]
N=3;
n=3;
X=[0;0;0]
gaussseidel(A,n,N,X,b)
5. Numerical Differentiation
a. Programming to obtain derivatives numerically.
x=[1.0 1.2 1.4 1.6 1.8 2.0 2.2];
y=[2.7182 3.3201 4.552 4.9530 6.0496 7.3891 9.0250];
c=1;
for i=1:6
d1(c)=y(i+1)-y(i);
c=c+1;
end
c=1;
for i=1:5
d2(c)=d1(i+1)-d1(i);
c=c+1;
end
c=1;
for i=1:4
d3(c)=d2(i+1)-d2(i);
c=c+1;
end
c=1;
for i=1:3
d4(c)=d3(i+1)-d3(i);
c=c+1;
end
c=1;
for i=1:2
d5(c)=d4(i+1)-d4(i);
c=c+1;
end
c=1;
for i=1:1
d6(c)=d5(i+1)-d5(i);
c=c+1;
end
x0=1.2
h=0.2;
f1=((d1(2)-d2(2)/2+d3(2)/3-d4(2)/4+d5(2)/5)/h);
printf('the first derivative of function at 1.2 is %f\n',f1)
f2=(d2(2)-d3(2)+(11*d4(2))/12-(5*d5(2))/6)/h^2;
printf('the second derivative of function at 1.2 is %f\n',f2)
6. Numerical Integration
a. Program for numerical integration using Trapezoidal rule.
function[i]=trapezoidal(a,b,n,f)
h=(b-a)/n;
x=(a:h:b);
y=f(x);
m=length(y);
i=y(1)+y(m);
for j=2:m-1;
i=i+2*y(j);
end;
i=(h/2)*i;
return(i);
endfunction
deff('[y]=f(x)','y=exp(x)')
trapezoidal(0,1,4,f)
b. Program for numerical integration using Simpson’s 1/3rd rule
function[i]=simpsons13(a,b,n,f)
h=(b-a)/n;
x=(a:h:b);
y=f(x);
m=length(y);
i=y(1)+y(m);
for j=2:m-1;
if(modulo(j,2)==0)then
i=i+4*y(j);
else
i=i+2*y(j);
end;
end;
i=(h/3)*i;
return(i);
endfunction
deff('[y]=f(x)','y=exp(x)')
simpsons13(0,4,4,f)
c. Program for numerical integration using Simpson’s 3/8th rule.
function[i]=simpsons38(a,b,n,f)
h=(b-a)/n;
x=(a:h:b);
y=f(x);
m=length(y);
i=y(1)+y(m);
for j=2:m-1
if(modulo(j-1,3)==0)then
i=i+2*y(j);
else
i=i+3*y(j);
end;
end;
i=(3*h/8)*i;
return(i);
endfunction
deff('[y]=f(x)','y=4+2*sin(x)')
simpsons38(0,%pi,6,f)
7. Solution of differential equations
a. Program to solve differential equation using Euler’s method
function [Y0]=eular(X0, Y0, h, yest, f)
n=(yest-X0)/h
for i=1:n
Y0=Y0+f(X0,Y0)*h;
X0=X0+h;
disp(Y0)
end;
endfunction
deff('[y]=f(a,b)','y=b-a*b+a');
eular(0,1,0.2,1,f)
b. Program to solve differential equation using modified Euler’s method.
function [y10]=eularmod(x0, y0, h, n, f)
x1=x0+h;
y10=y0+h*f(x0,y0)
while(n>1)
x0=x0+h;
x1=y10;
y10=y0+(h/2)*(f(x0,y0)+f(x1,y10));
if(abs(y10-x1)<0.001)
y10
abort;
end
n=n-1;
y10
end
endfunction
deff('[y]=f(a,b)','y=log(a+b)');
eularmod(1,2,0.2,10,f)
c. Program to solve differential equation using Runge-kutta 2nd order and 4th order methods.
2nd order
function [y]=f(a, b)
y=b-a;
endfunction
x0=0;
y0=2;
h=0.1;
for n=1:4
k1=h*f(x0,y0);
k2=h*f(x0+h,y0+k1);
y0=y0+(k1+k2)/2;
x0=x0+h;
printf('values of x0=%g\t and y0=%g\n',x0,y0);
end
4th order
function [y]=f(a, b)
y=b-a;
endfunction
x0=0;
y0=1;
h=0.25;
for n=1:4
k1=h*f(x0,y0);
k2=h*f(x0+h/2,y0+k1/2);
k3=h*f(x0+h/2,y0+k2/2);
k4=h*f(x0+h,y0+k3);
y0=y0+(k1+2*k2+2*k3+k4)/6;
x0=x0+h;
printf('values of x0=%g\t and y0=%g\n',x0,y0);
end
8. Regression
a. Program for Linear regression.
x=[1,2,3,4,5,6,7];
y=[0.5,2.5,2,4,3.5,6,5.5];
n=7;
s=0;
xsq=0;
xsum=0;
ysum=0;
for i =1.7
s=s+(det(x(1,i)))*(det(y(1,i)));
xsq=xsq+(det(x(1,i))^2);
xsum=xsum+det(x(1,i));
ysum=ysum+det(y(1,i));
end
disp(s,"sum of product of x and y=")
disp(xsq,"sum of square of x=")
disp(xsum,"sum of all the x=")
disp(ysum,"sum of all the y=")
a=xsum/n;
b=xsum/n;
a1=(n*s-xsum*ysum)/(n*xsq-xsum^2);
a0=b-a*a1;
disp(a1,"a1=")
disp(a0,"a0=")
disp("The equation of the line obtained is y=a0+a1*x")
b. Program for Polynomial Regression.
x=[0,1,2,3,4,5];
y=[2.1,7.7,13.6,27.2,40.9,61.1];
sumy=0;
sumx=0;
m=2;
n=6;
s=0;
xsqsum=0;
xcsum=0;
x4sum=0;
xysum=0;
x2ysum=0;
rsum=0;
usum=0;
for i=1:6
sumy=sumy+y(i);
sumx=sumx+x(i);
r(i)=(y(i)-s/n)^2;
xsqsum=xsqsum+x(i)^2;
xcsum=xcsum+x(i)^3;
x4sum=x4sum+x(i)^4;
xysum=xysum+x(i)*y(i);
x2ysum=x2ysum+y(i)*x(i)^2;
rsum=r(i)+rsum;
end
disp(sumy,"sum y=")
disp(sumx,"sum x")
xavg=sumx/n;
yavg=sumy/n;
disp(xavg,"xavg=")
disp(yavg,"yavg=")
disp(xsqsum,"sum x^2=")
disp(xcsum,"sum x^3=")
disp(x4sum,"sum x^4=")
disp(xysum,"sum x*y=")
disp(x2ysum,"sum x^2*y=")
J=[n,sumx,xsqsum;sumx,xsqsum,xcsum;xsqsum,xcsum,x4sum];
I=[sumy;xysum;x2ysum];
X=inv(J)*I;
a0=det(X(1,1));
a1=det(X(2,1));
a2=det(X(3,1));
for i=1:6
u(i)=(y(i)-a0-a1*x(i)-a2*x(i)^2)^2;
usum=usum+u(i);
end
disp(r,"(yi-yavg)^2=")
disp(u,"(yi-a0-a1*x-a2*x^2)^2=")
plot(x,y);
xtitle('x vs y','x','y');
syx=(usum/(n-3))^0.5;
disp(syx,"The standard error of the estimate based on regression polynomial=")
R2=(rsum-usum)/(rsum);
disp("%",R2*100,"Percentage of original uncertainty that has been explained by the model=")
c. Program for multiple linear regression.
x1=[0,2.2,2.5,1,4,7];
x2=[0,1,2,3,6,2];
x1sum=0;
x2sum=0;
ysum=0;
x12sum=0;
x22sum=0;
x1ysum=0;
x2ysum=0;
x1x2sum=0;
n=6;
for i=1:6
y(i)=5+4*x1(i)-3*x2(i);
x12(i)=x1(i)^2;
x22(i)=x2(i)^2;
x1x2(i)=x1(i)*x2(i);
x1y(i)=x1(i)*y(i);
x2y(i)=x2(i)*y(i);
x1sum=x1sum+x1(i);
x2sum=x2sum+x2(i);
ysum=ysum+y(i);
x1ysum=x1ysum+x1y(i);
x2ysum=x2ysum+x2y(i);
x1x2sum=x1x2sum+x1x2(i);
x12sum=x12sum+x12(i);
x22sum=x22sum+x22(i);
end
X=[n,x1sum,x2sum;x1sum,x12sum,x1x2sum;x2sum,x1x2sum,x22sum];
Y=[ysum;x1ysum;x2ysum];
Z=inv(X)*Y;
a0=det(Z(1,1));
a1=det(Z(2,1));
a2=det(Z(3,1));
disp(a0,"a0=")
disp(a1,"a1=")
disp(a2,"a2=")
disp("Thus,y=a0+a1*x1+a2*x2")
d. Program for non-linear regression.
x=[0.25,0.75,1.25,1.75,2.25];
y=[0.25,0.57,0.68,0.74,0.79];
a0=1
a1=1;
sr=0.0248;
for i=1:5
pda0(i)=1-exp(-a1*x(i));
pda1(i)=a0*x(i)*exp(-a1*x(i));
end
Z0=[pda0(1),pda1(1);pda0(2),pda1(2);pda0(3),pda1(3);pda0(4),pda1(4);pda0(5),pda1(5)]
disp(Z0,"Z0=")
R=Z0'*Z0;
S=inv(R);
for i=1:5
y1(i)=a0*(1-exp(a1*x(i)));
D(i)=y(i)-y1(i);
end
disp(D,"D=")
M=Z0'*D;
X=S*M;
disp(X,"X=")
a0=a0+det(X(1,1));
a1=a1+det(X(2,1));
disp(a0,"The value of a0 after 1st iteration=")
disp(a1,"The value of a1 after 1st iteration")
9. Random variables and distributions
a. Program to generate random variables.
a=3;
x=27;
c=2;
m=10;
for i=0:m
xi=modulo((a*x+c),m)
printf("%f,",xi/m)
x=xi
end
b. Program to fit binomial distribution.
function [expfre]=binfit(x, obsfre)
n=length(x)-1;
N=sum(obsfre);
xbar=sum(x*obsfre')/N;
p=xbar/n;
q=1-p;
for r=0:n
prob(r+1)=factorial(n)*p^r*q^(n-r)/(factorial(r)*factorial(n-r));
end
expfre=round(prob*N);
printf('expected frequencies\n');
printf('------------------------------\n');
return(expfre);
endfunction
x=[0,1,2,3,4,5];
obsfre=[12,25,36,21,10,8];
binfit(x,obsfre)
c. Program to fit Poisson distribution.
function [expfre]=poissonfit(x, obsfre)
n=length(x)-1;
N=sum(obsfre);
xbar=x*obsfre'/N;
m=xbar;
prob(1)=exp(-m);
for i=1:n
prob(i+1)=m*prob(i)/(i)
end
expfre=round(prob*N);
printf('expected frequenics\n');
printf('------------------------------\n');
return(expfre);
endfunction
x=[0,1,2,3,4,5];
obsfre=[100,70,45,20,10,5];
poissonfit(x,obsfre)
10. Distributions
a. Program for Uniform distribution.
function [expfre]=unifit(x, obsfre)
n=length(x)-1;
N=sum(obsfre);
xbar=sum(x*obsfre')/N;
prob(1)=1/(n+1);
for r=0:n
prob(r+1)=1/(n+1);
end
expfre=round(prob*N);
printf('expected frequencies\n');
printf('------------------------------\n');
return(expfre);
endfunction
x=[0,1,2,3,4,5];
obsfre=[12,50,100,80,45,5];
unifit(x,obsfre)
b. Program for Bernoulli distribution
function [expfre]=bernfit(x, obsfre)
n=length(x)-1;
N=sum(obsfre);
xbar=sum(x*obsfre')/N;
p=xbar/n;
q=1-p;
for r=0:n
prob(r+1)=p^r*q^(n-r);
end
expfre=round(prob*N);
printf('expected frequencies\n');
printf('------------------------------\n');
return(expfre);
endfunction
x=[0,1];
obsfre=[12,25];
bernfit(x,obsfre)
c. Program for Negative binomial distribution.
function [expfre]=negbinfit(x, obsfre)
n=length(x)-1;
N=sum(obsfre);
xbar=sum(x*obsfre')/N;
rawmom2=sum(obsfre*(x^2)')/N;
var=rawmom2-xbar^2;
p=xbar/var;
q=1-p;
r=round(p*xbar/q);
prob(1)=p^r;
for i=1:n
prob(i+1)=((i+1+r)/(i-1+1))*q*prob(i)
end
printf('expected frequencies\n');
printf('------------------------------\n');
expfre=round(prob*N);
return(expfre);
endfunction
x=[0,1,2,3,4,5];
obsfre=[12,50,100,80,45,5];
negbinfit(x,obsfre)