MATLAB Program No: 01
Program Name: An M-File to implement the Bisection Method.
Problem: Find a real root of the equation f(x)=x3-x-1=0 by using bisection method.
function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)
if nargin<3,error('at least 3 input arguments required'),end
test = func(xl,varargin{:})*func(xu,varargin{:});
if test>0,error('no sign change'),end
if nargin<4|isempty(es), es=0.0001;end
if nargin<5|isempty(maxit), maxit=50;end
iter = 0; xr = xl; ea = 100;
while (1)
xrold = xr;
xr = (xl + xu)/2;
iter = iter + 1;
if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end
test = func(xl,varargin{:})*func(xr,varargin{:});
if test < 0
xu = xr;
elseif test > 0
xl = xr;
else
ea = 0;
end
if ea <= es | iter >= maxit,break,end
end
root = xr; fx = func(xr, varargin{:});
Input program: y=@(x)x^3-2*x-5;
dy=@(x)3*x^2-2;
[root ea iter]=newtraph(y,dy,2,0.00001)
OUTPUT
On Command Window On Workspace
MATLAB Program No: 02
Program Name: An M-File to implement the Newton Raphson Method.
Problem: Use the Newton-Raphson method to find a root of the equation x3-2x-5=0
function
[root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,varargin)
if nargin<3,error('at least 3 input arguments required'),end
if nargin<4|isempty(es),es=0.0001;end
if nargin<5|isempty(maxit),maxit=50;end
iter = 0;
while (1)
xrold = xr;
xr = xr - func(xr)/dfunc(xr);
iter = iter + 1;
if xr ~= 0, ea = abs((xr - xrold)/xr) * 100; end
if ea <= es | iter >= maxit, break, end
end
root = xr;
Input program: y=@(x)x^3-2*x-5;
dy=@(x)3*x^2-2;
[root ea iter]=newtraph(y,dy,2,0.00001)
OUTPUT
On Command Window On Workspace
MATLAB Program No: 03
Program Name: An M-File to implement Linear regression.
Problem: An experiment gave the following table of values for the dependent variable y for
a set of known values of x. Obtain an appropriate least squares fit for the data and plot y vs x
with MATLAB.
x 1 2 3 4 5 6 7 8 9
y 5.5 7.0 9.6 11.5 12.6 14.4 17.6 19.5 20.5
function [a, r2] = linregr(x,y)
n = length(x);
if length(y)~=n, error('x and y must be same length'); end
x = x(:); y = y(:); % convert to column vectors
sx = sum(x); sy = sum(y);
sx2 = sum(x.*x); sxy = sum(x.*y); sy2 = sum(y.*y);
a(1) = (n*sxy-sx*sy)/(n*sx2-sx^2);
a(2) = sy/n-a(1)*sx/n;
r2 = ((n*sxy-sx*sy)/sqrt(n*sx2-sx^2)/sqrt(n*sy2-sy^2))^2;
% create plot of data and best fit line
xp = linspace(min(x),max(x),2);
yp = a(1)*xp+a(2);
plot(x,y,'o',xp,yp)
grid on
Input program: x=[1 2 3 4 5 6 7 8 9];
y=[5.5 7 9.6 11.5 12.6 14.4 17.6 19.5 20.5];
linregr(x,y)
OUTPUT
On Command Window Plot y vs x;
MATLAB Program No: 04
Program Name: An M-File to implement Lagrange Interpolation.
Problem: Given the following tables of values
x 0.4 0.5 0.7 0.8
y -0.916 -0.693 -0.357 -0.223
Use the above m-file to find the value of f(0.6).
function yint = Lagrange(x,y,xx)
n = length(x);
if length(y)~=n, error('x and y must be same length'); end
s = 0;
for i = 1:n
product = y(i);
for j = 1:n
if i ~= j
product = product*(xx-x(j))/(x(i)-x(j));
end
end
s = s+product;
end
yint = s;
Input program: x=[.4 .5 .7 .8];
y=[-.916 -.693 -.357 -.223];
fx=Lagrange(x,y,0.6)
OUTPUT
On Command Window On Workspace
MATLAB Program No: 05
Program Name: An M-File to implement Newton interpolation.
Problem: Values of x (in degrees) and sin x are given in the following table:
x (in degrees) 15 20 25 30 35 40
sin x 0.258819 0.3420201 0.4226183 0.5 0.5735764 0.6427876
Determine the value of sin 38°.
function yint = Newtint(x,y,xx)
n = length(x);
if length(y)~=n, error('x and y must be same length'); end
b = zeros(n,n);
% assign dependent variables to the first column of b.
b(:,1) = y(:); % the (:) ensures that y is a column vector.
for j = 2:n
for i = 1:n-j+1
b(i,j) = (b(i+1,j-1)-b(i,j-1))/(x(i+j-1)-x(i));
end
end
% use the finite divided differences to interpolate
xt = 1;
yint = b(1,1);
for j = 1:n-1
xt = xt*(xx-x(j));
yint = yint+b(1,j+1)*xt;
end
Input program: x=[15 20 25 30 35 40];
y=[.2588190 .3420201 .4226183 .5 .5735764 .6427876];
Newtint(x,y,38)
OUTPUT
On Command Window On Workspace
MATLAB Program No:06
Program Name: An M-File to implement the Trapezoidal Rule.
1
1
Problem: Evaluate I= dx correct to three decimal places with h=0.125.
0 1 x
function I = trap(func,a,b,n,varargin)
if nargin<3,error('at least 3 input arguments required'),end
if ~(b>a),error('upper bound must be greater than lower'),end
if nargin<4|isempty(n),n=1000000;end
x = a; h = (b - a)/n;
s=func(a,varargin{:});
for i = 1 : n-1
x = x + h;
s = s + 2*func(x,varargin{:});
end
s = s + func(b,varargin{:});
I = (b - a) * s/(2*n);
Input program: y=@(x)1/(1+x);
trap(y,0,1,8)
OUTPUT
On Command Window On Workspace
MATLAB Program No: 07
Program Name: An M-File to implement the Trapezoidal Rule for unequally spaced
data.
function I = trapuneq(x,y)
% trapuneq: unequal spaced trapezoidal rule quadrature
% I = trapuneq(x,y):
% Applies the trapezoidal rule to determine the integral
% for n data points (x, y) where x and y must be of the
% same length and x must be monotonically ascending
% input:
% x = vector of independent variables
% y = vector of dependent variables
% output:
% I = integral estimate
if nargin<2,error('at least 2 input arguments required'),end
if any(diff(x)<0),error('x not monotonically ascending'),end
n = length(x);
if length(y)~=n,error('x and y must be same length'); end
s = 0;
for k = 1:n-1
s = s + (x(k+l)-x(k))*(y(k)+y(k+l))/2;
end
I = s;
Input program:
x = [0 .12 .22 .32 .36 .4 .44 .54 .64 .7 .8];
y = 0.2+25*x-200*x.^2+675*x.^3-900*x.^4+400*x.^5;
trapz(x,y)
OUTPUT
On Command Window On Workspace
