MATLAB Programs on Unit 3 & 4
Program no 6
MATLAB Program for Regula-Falsi/Newton-Raphson Method
‘Regula-Falsi Method’
%Solution of Algebraic and Transcendental Equation by Regula-Falsi Method
clc;
close all;
clear all;
syms x;
f=input('Enter the function:\n'); %Enter the Function here
n=input('Enter the number of decimal places:\n');
epsilon = 10^-(n+1); %Error tolerance
x0 = input('Enter the lower limit:\n');
x1 = input('Enter the upper limit:\n');
for i=1:100
f0=vpa(subs(f,x,x0)); %Calculating the value of function at x0
f1=vpa(subs(f,x,x1)); %Calculating the value of function at x1
y=x1-((x1-x0)/(f1-f0))*f1; %[x0,x1] is the interval of the root
err=abs(y-x1);
if err<epsilon %checking the amount of error at each iteration
break
end
f2=vpa(subs(f,x,y));
if (f1)*(f2)<1
x0=y; %taking the next interval as[x0,x1] = [y,x1]
x1=x1;
else
x0=x0; %taking the next interval as[x0,x1] = [x0,y]
x1=y;
end
end
y = y - rem(y,10^-n); %Displaying upto required decimal places
fprintf('The Root is : %f \n',y);
‘Newton-Raphson Method’
%Solution of Algebraic and Transcendental Equation by Newton-Raphson Method
clc;
close all;
clear all;
syms x;
f=input('Enter the function:\n'); %Enter the Function here
g=diff(f); %The Derivative of the Function
n=input('Enter the number of decimal places:');
epsilon = 5*10^-(n+1)
x0 = input('Enter the intial approximation:');
for i=1:100
f0=vpa(subs(f,x,x0)); %Calculating the value of function at x0
f0_der=vpa(subs(g,x,x0)); %Calculating the value of function derivative at x0
y=x0-f0/f0_der; % The Formula
err=abs(y-x0);
if err<epsilon %checking the amount of error at each iteration
break
end
x0=y;
end
y = y - rem(y,10^-n); %Displaying upto required decimal places
fprintf('The Root is : %f \n',y);
fprintf('No. of Iterations : %d\n',i);
Program no 7
Interpolation using Newton’s Forward/Backward Difference Formula
‘Forward Difference Formula’
%Newtonns Forward Interpolation Formula MATLAB Program
clc;
close all;
clear all;
x = input('Enter the values of x:\n'); % inputting values of x
fx = input('Enter the values of y:\n'); % inputting values of y
dt=zeros(10,10); % function
for i=1:6 dt(i,1)=x(i); % for loop
dt(i,2)=fx(i); % calling function
end
n=5; % number of iterations
for j=3:10
for i=1:n
dt(i,j)=dt(i+1,j-1)-dt(i,j-1)
end
n=n-1;
end
h=x(2)-x(1) % finding the value of h
xp=1.5; % defining the value of xp
for i=1:5
q=(xp-x(i))/h; % calculating number of intervals
if (q>0 && q<1)
p=q;
end
end
p
l=xp-(p*h)
for i=1:5
if(l==x(i))
r=i;
end
end % calculating different value of y
f0=fx(r);
f01=dt(r,3);
f02=dt(r,(3+1));
f03=dt((r),(3+2));
f04=dt((r),(3+3));
% using the forward interpolation formula
fp=(f0)+((p*f01)+(p*(p-1)*f02)/(2)) + ((p*(p-1)*(p-2)*f03)/(6))+((p*(p-1)*(p-
2)*(p-3)*f04)/(24))
‘Backward Difference Formula’
%Newtons Backward Interpolation Formula MATLAB Program
clc;
close all;
clear all;
x = input('Enter the values of x:\n'); % inputting the values of x
fx = input('Enter the values of y:\n'); % inputting the value of y
dt=zeros(7,7); % declaring function
for i=1:6 % stating loop
dt(i,1)=x(i);
dt(i,2)=fx(i);
end
n=5;
for j=3:7
for i=1:n % using for loop
dt(i,j)=dt(i+1,j-1)-dt(i,j-1) % defining dt
end
n=n-1;
end
h=x(2)-x(1) % finding the value of h
xp=27; % defining xp
for i=1:6
q=(xp-x(i))/h;
if (q>0&&q<1)
p=q;
end
end
p
l=xp-(p*h)
for i=1:6
if(l==x(i))
r=i+1;
end
end
% finding different value of y
f0=fx(r);
f01=dt((r-1),3);
f02=dt((r-2),(3+1));
f03=dt((r-3),(3+2));
f04=dt((r-4),(3+3));
%using backward difference formula
fp=(f0)+((p*f01)+(p*(p+1)*f02)/(2)) + ((p*(p+1)*(p+2)*f03)/(6))
Program no 8
Computation of Area under the curve using Trapezoidal, Simpson’s
1/3rd and 3/8th Rule
%Numerical Integration
clc;
close all;
clear all;
syms x f(x)
f(x) = input('Enter the integrand:\n');
x0 = input('Enter the lower limit:\n');
xn = input('Enter the upper limit:\n');
n = input('Enter the number of intervals:\n');
h = (xn-x0)/n;
xinit = x0;
for i=1:n+1
x(i) = xinit;
xinit = x(i) + h;
end
for i=1:n+1
y(i) = subs(f,x(i));
end
%Trapezoidal Rule
Area_by_Trapezoidal_Rule = (h/2)*((y(1)+y(n+1))+2*(sum(y)-y(1)-y(n+1)))
%Simpson's 1/3rd Rule
odd = 0;
even = 0;
for i=2:2:n
odd = odd + y(i);
end
for i = 3:2:n
even = even + y(i);
end
Area_by_1_3rd_Rule = (h/3)*((y(1)+y(n+1))+4*(odd)+2*(even))
%Simpson's 3/8th Rule
nonthird = 0;
third = 0;
for i = 1:1:n-1
if rem(i,3)==0
third = third + y(i+1);
else
nonthird = sum(y)-y(1)-y(n+1)-third;
end
end
Area_by_3_8th_Rule = (3*h/8)*((y(1)+y(n+1))+3*nonthird+2*third)
Program no 9
Solution of ODE of first order and first degree by Taylor’s series and
Modified Euler’s Method
‘Taylor’s Series Method’
%Taylor's Series Method for f(x,y)=x^2*y-1
clear
clc
x0 = input('Enter initial value of x:\n');
y0 = input('Enter initial value of y:\n');
xn = input('Enter the value of x at which y is to be calculated:\n');
h = input('Enter the step size:\n');
n = (xn-x0)/h;
for i = 1:n
ysol = y0 + h*f(x0,y0) + ((h^2)/2)*f1(x0,y0) + ((h^3)/6)*f2(x0,y0); %Taylor's
series upto 2nd order
y0 = ysol;
x0 = x0 + h;
end
fprintf('The solution after %f iterations is: %f',n,ysol)
function dy = f(x,y)
dy = x^2*y-1; %Enter your f(x,y)
end
function d2y = f1(x,y)
d2y = (2*x+x^4)*y-x^2; %Enter your f'(x,y)
end
function d3y = f2(x,y)
d3y = x^6*y-x^4+6*x^3*y-4*x+2*y; %Enter your f''(x,y)
end
‘Modified Euler’s Method’
%Modified Euler's Method written for f(x,y) = x+y.
clear
clc
x0 = input('Enter initial value of x : \n');
y0 = input('Enter initial value of y : \n');
x1 = input('Enter the value of x at which y is to be calculated : \n');
h = input('Enter the value of h : \n');
n = (x1 - x0)/h;
for i = 1:n
y1p = y0 + h*f(x0,y0);
y1c = y0 + (h/2)*(f(x0,y0)+f(x0+h,y1p));
err = abs(y1p - y1c);
while err > 0.00001
y1c = y0 + (h/2)*(f(x0,y0)+f(x0+h,y1p));
err = abs(y1p - y1c);
y1p = y1c;
end
y0 = y1c;
x0 = x0 + h;
end
fprintf('The solution after %f iterations is: %f',n,y1c)
function z = f(x,y)
z = x+y; %Enter your f(x,y)
end
Program no 10
Solution of first order and first degree ODE by Runge-Kutta 4th order
and Milne’s method
‘Runge-Kutta 4th Order Method’
%Runge-Kutta 4th Order Method for f(x,y)=x+y
clear
clc
x0 = input('Enter the initial value of x:\n');
y0 = input('Enter the initial value of y:\n');
xn = input('Enter the value of x at which y is to be calculated:\n');
h = input('Enter the step size:\n');
n = (xn - x0)/h;
%Runge-Kutta Formula
for i =1:n
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);
yn = y0 + (1/6)*(k1+(2*k2)+(2*k3)+k4);
x0 = x0 + h;
y0 = yn;
end
fprintf('The solution after %f iterations is: %f',n,yn)
function z = f(x,y)
z = x+y;
end
‘Milne’s Method’
%Milne's Method for f(x,y)=x-y^2
clear
clc
x0 = input('Enter initial value of x:\n');
y0 = input('Enter initial value of y:\n');
xn = input('Enter the value of x at which y is to be calculated:\n');
n = input('Enter the number of steps:\n');
h = (xn-x0)/n;
%Taylor's series to get the ordinates (y(i))
for i = 1:n
ysol(i) = y0(i) + h*f(x0(i),y0(i)) + ((h^2)/2)*f1(x0(i),y0(i)) +
((h^3)/6)*f2(x0(i),y0(i));
y0(i+1) = ysol(i);
x0(i+1) = x0(i) + h;
end
%Milne's Method from here
for i = 1:n
F(i) = f(x0(i),y0(i));
end
for i = 5:n
y0(i) = y0(i-4) + (4*h/3)*(2*F(i-3)-F(i-2)+2*F(i-1)); %Predictor Formula
F(i) = f(x0(i),y0(i));
y0(i) = y0(i-2) + (h/3)*(F(i-2)+4*F(i-1)+F(i)); %Corrector Formula
end
fprintf('The solution is: %f\n',y0(n+1))
function dy = f(x,y)
dy = x-y^2; %Enter f(x,y)
end
function d2y = f1(x,y)
d2y = 1-2*y*(x-y^2); %Enter f'(x,y)
end
function d3y = f2(x,y)
d3y = -6*y^4+8*x*y^2-2*y-2*x^2; %Enter f''(x,y)
end
