%PROGRAM FOR BISECTION METHOD clc; % clears the command window clear all all; ; % clears the workspace f=input('give f=input('give f(x)= '); '); % syntax for input n=input('give n=input('give no. of iterations, n= '); '); c=1; % initialisation of while loop while (c while (c >= 0) a=input('give a=input('give a= '); '); b=input('give b=input('give b= '); '); c=f(a)*f(b); end ; for i=1:n for i=1:n x(i)=(a+b)/2; g=f(x(i)); if (f(a)*g if (f(a)*g < 0) b=x(i); else a=x(i); end ; disp(x(i)); end ;
OUTPUT: give f(x)= @(x)(x^3-x-1) give no. of iteations, n= 8 give a= 1 give b= 2 1.5000 1.2500 1.3750 1.3125 1.3438 1.3281 1.3203 1.3242
SOLVER: >> p=[1 0 -1 -1]; >> r=roots(p) r = 1.3247 -0.6624 + 0.5623i -0.6624 - 0.5623i
% program for false position method clc; % clears the command window clear all all; ; % clears the workspace f=input('give f=input('give f(x)= '); '); % syntax for input n=input('give n=input('give no. of iterations, n= '); '); c=1; % initialisation of while loop while (c while (c >= 0) a=input('give a=input('give a= '); '); b=input('give b=input('give b= '); '); c=f(a)*f(b); end ; for i=1:n for i=1:n x(i)=(a*f(b)-b*f(a))/(f(b)-f(a)); if (f(a)*f(x(i))< if (f(a)*f(x(i))< 0) b=x(i); else a=x(i); end ; disp(x(i)); end ;
Output: give f(x)= @(x)(x^3-5*x+3) give no. of iterations, n= 6 give a= 0 give b= 1 0.7500 0.6761 0.6604 0.6573 0.6568 0.6566
Solver: >> p=[1 0 -5 3]; >> r=roots (p)
% program for successive approximation method clc; clear all; f=input('give f(x)= '); n=input('no. of iterations n= '); c=1; % for initialising the while loop while (c >= 0) a=input('give intial assumption a= '); b=input('give final assumption b= '); c=f(a)*f(b); end ; x0=(a+b)/2; e=2; while (e >=1) phi=input('give phi(x)= '); phid=input('give phi dash(x)= '); e = abs(phid(x0)); end ; for i=1:n x1=phi(x0); disp(x1); x0=x1; end ; Output: give f(x)= @(x)(x^3-x-1) no. of iterations n= 5 give intial assumption a= 0 give final assumption b= 1 give intial assumption a= 1 give final assumption b= 2 give phi(x)= @(x)(x^3-1) give phi dash(x)= @(x)(3*x^2) give phi(x)= @(x)((1+x)/x^2) give phi dash(x)= @(x)((-2/x^3)-(1/x^2)) give phi(x)= @(x)((1+x)^(1/3)) give phi dash(x)= @(x)((1/3)*(1+x)^(-2/3)) 1.3572 1.3309 1.3259 1.3249 1.3248 Solver: >> p=[1 0 -1 -1]; >> r=roots(p) r = 1.3247 -0.6624 + 0.5623i -0.6624 - 0.5623i
% Newton Raphson Method to find roots of given equation clc; clearall; f=input ('give f(x)= '); fd=input('give fdash(x)= '); n=input('give no. of iterations n= '); c=1;
% initialisation to enter in while loop for first t
while (c > 0) a=input('give initial assumption a= '); b=input('give final assumption b= '); c=f(a)*f(b); end ; x0=(a+b)/2; for i=1:n x1=x0-(f(x0)/fd(x0)); x0=x1; disp(x1); end ;
Output: give f(x)= @(x)(exp(x)*cos(x)-1.4) givefdash(x)= @(x)(exp(x)*(cos(x)-sin(x))) give no. of iterations n= 3 give initial assumption a= 0 give final assumption b= 1 0.4286 0.4335 0.4336
Solver: >>syms x >> p=(exp(x)*cos(x)-1.4); >>root=solve(p) root = .43356087535265770522790906547161
% Trapezoidal integration clc; clearall; f=input('give f(x)= '); x0=input('give initial limit,x0= '); xn=input('give final limit,xn= '); n=input('give steps,n= '); h=(xn-x0)/n; for i=0:n x(i+1)=x0+(i)*h; y(i+1)=f(x(i+1)); end ; a=y(1)+y(n+1); b=0; for i=2:n b=b+y(i); end ; I=(h/2)*(a+2*b); fprintf('integration I= %f \n',I)
Output: give f(x)= @(x)(4*x+2) give initial limit,x0= 3 give final limit,xn= 5 givesteps,n= 4 integration I= 36.000000
Solver: >> x=(3:0.5:5); >> I=trapz(x,4*x+2) I = 36
%Program for Simpson’s 1/3 Rule clc; clearall; f=input('Give f(x) = '); x0=input('Give first limit, x0 = '); xn=input('Give final limit, xn = '); n=input('Give no. of steps that divisible by 2, n = '); h=(xn-x0)/n; ans=0; for i=1:n-1 if (mod(i,2)==0) ans=ans+2*f(x0+(i*h)); else ans=ans+4*f(x0+(i*h)); end ; end ; ans1=f(x0)+f(xn); area=ans+ans1; area=area*(h/3); disp (area);
Output: Give f(x) = @(x)(log(x+1)+sin(2*x)) Give first limit, x0 = 0 Give final limit, xn = 0.8 Give no of steps that divisible by 2, n = 8 0.7726
Solver: >> f='(log(x+1)+sin(2*x))'; >> I=quad(f,0,0.8) I = 0.7726
%Program for Simpson’s 3/8 Rule clc; clearall; f=input('Give f(x) = '); x0=input('Give lower limit, x0 = '); xn=input('Give upper limit, xn = '); n=input('Give no. of steps that is divisible by 3, n = '); h=(xn-x0)/n; for i=1:n+1 x(i)=x0+(i-1)*h; y(i)=f(x(i)); end ; ans=y(1)+y(n+1); for i=1:n-1 if (mod(i,3)==0) ans=ans
+ 2*y(i+1);
else ans=ans + 3*y(i+1); end ; end ; area=ans*(3*h/8); disp (area);
Output: Give f(x) = @(x)(4*x-1) Give lower limit, x0 = 1 Give upper limit, xn = 4 Give no. of steps that is divisible by 3, n = 6 27
Solver: >> f='4*x-1'; >> I=quad (f,1,4) I = 27
% Program For Lagranges Interpolation clc; clearall; n=input('enter total no. of data points, n= '); x=input('enter x= '); y=input('enter y= '); xg=input('enter x for which y is unknown, xg= '); yg=0; for j=1:n p=1; q=1; for k=1:n if j~=k p=p*(xg-x(k)); q=q*(x(j)-x(k)); end ; end ; L(j)=p*y(j)/q; yg=yg+L(j); end ; fprintf('at x=%f, y=%f',xg,yg);
Output: enter total no. of data points, n= 5 enter x= [1 2 3 4 5] enter y= [20 40 60 80 100] enter x for which y is unknown, xg= 2.5 at x=2.500000, y=50.000000
Solver: >> x=[ 1 2 3 4 5]; >> y=[20 40 60 80 100]; >>yg= interp1(x,y,2.5) yg = 50
% Newton forward difference interpolation x=input('enter all values for x='); y=input('enter all values for y='); n=input('enter total no. of data points n='); xi=input('enter x for which y is unknown xi='); h=x(2)-x(1); u=(xi-x(1))/h; for i=2:n for j=1:(n-(i-1)) y(i,j)=y(i-1,j+1)-y(i-1,j); end end r=y(1,1); for i=2:n d=y(i,1); for k=0:(i-2) d=d*(u-k); end d=d/factorial(i-1); r=r+d; end disp(r);
output: enter all values for x=[1 2 3 4 5] enter all values for y=[10 20 30 40 50] enter total no. of data points n=5 enter x for which y is unknown xi=1.5 15
Solver: >> x= [1 2 3 4 5]; >> y= [10 20 30 40 50]; >> Y=interp1(x,y,1.5) Y =
15
% Gauss Seidal Method to solve set of simultaneous equations clc; clearall; n=input('enter no.of unknowns, n= '); a=input('enter the set of eqn. in augmd. matrix form,a=\n ') for i=1:n x(i)=0; end ; for k=1:5 for i=1:n p=0; for j=1:n if i~=j p=p+a(i,j)*x(j); end ; end ; x(i)=(a(i,n+1)-p)/a(i,i); fprintf('x(%d,%d)=%f \n',k,i,x(i)); end ; end ;
Output: enter no.of unknowns, n= 3 enter the set of eqn. in augmd.matrix form,a= [4 1 1 5;1 6 2 19;-1 -2 -5 10;] x(1,1)=1.250000 x(1,2)=2.958333 x(1,3)=-3.433333 x(2,1)=1.368750 x(2,2)=4.082986 x(2,3)=-3.906944 x(3,1)=1.205990 x(3,2)=4.267983 x(3,3)=-3.948391 x(4,1)=1.170102
Solver: >> a=[4 1 1;1 6 2;-1 -2 -5]; >> b=[5; 19; 10]; >>linsolve(a,b) ans = 1.1649 4.2887 -3.9485
% Program On straight Line equation clc; clearall; n=input('enter total no. of data points ,n= '); x=input('enter set of x= '); y=input('enter set of y= '); sx=0; sy=0; ssq=0; sxy=0; for i=1:n sq(i)=x(i)*x(i); xy(i)=x(i)*y(i); sx=sx + x(i); sy=sy + y(i); ssq=ssq + sq(i); sxy=sxy + xy(i); end ; c=[n sxsy; sxssqsxy]; for j=1:3 c(2,j)=c(2,j)-(sx/n)*c(1,j); end ; b=c(2,3)/c(2,2); a=(sy-b*sx)/n; fprintf('the reqd. equ of line that fits best to given poin %f + %f x \n' ,a,b);
Output: enter total no. of data points ,n= 7 enter set of x= [1 2 3 4 5 6 7] enter set of y= [0.5 2.5 2 4 2.5 6 5.5] the reqd. equ of line that fits best to given points is y= 0.071429 + 0.803571 x
% program for least square approx. for fitting a second deg equation clc; clear all; close all; n= input('total no. of data points, n= '); x= input('enter array of x= '); y= input('enter array of y= '); sx=0; sy=0; sx2=0; sx3=0; sx4=0; sxy=0; sx2y=0; for i=1:n sx=sx + x(i); sy=sy + y(i); sx2=sx2 + x(i)^2; sx3=sx3 + x(i)^3; sx4=sx4 + x(i)^4; sxy=sxy + x(i)*y(i); sx2y=sx2y + (x(i)^2)*y(i); end ; m=[n sx sx2 sy; sx sx2 sx3 sxy; sx2 sx3 sx4 sx2y]; for i=1:3 a(i)=0; end ; for k=1:15 for i=1:3 p=0; for j=1:3 if(i~=j) p=p+m(i,j)*a(j); end ; end ; a(i)=(m(i,4)-p)/m(i,i); end ; end ; fprintf ('equation of given curve is \n y=%f + %f(x) + %f(x^2)\n',a(1),a(2),a(3));
Output:
Solver: >> x=[-3 -2 -1 0 1 2 3]; >> y=[12 4 1 2 7 15 30]; >>polyfit(x,y,2) ans = 2.1190
2.9286
1.6667
% Program on Euler’s method clc; clearall; f=input('enter given function f= '); x(1)=input('enter the x0= '); xn=input('enter the xn= '); y(1)=input('enter the y0= '); h=input('enter the increment in x, h= '); n=(xn-x(1))/h; for i=1:n+1 x(i+1)=x(i)+h; y(i+1)=y(i)+h*f(x(i),y(i)); fprintf('\n at x(%d)=%f value of y(%d)=%f',i+1,x(i+1),i+1,y end ;
Output: enter given function f= @(x,y)(x+2*y) enter the x0= 1 enter the xn= 1.4 enter the y0= 1 enter the increment in x, h= 0.1 at x(2)=1.100000 value of y(2)=1.300000 at x(3)=1.200000 value of y(3)=1.670000 at x(4)=1.300000 value of y(4)=2.124000 at x(5)=1.400000 value of y(5)=2.678800
Solver: Save .m file with name EULRsolver functiondy=f(x,y) dy=(x+2*y); @ command window: >> [x,y]=ode45('EULRsolver',[1:0.1:1.4],1) x =
1.0000 1.1000 1.2000 1.3000
% PROGRAM FOR MODIFIED EULERS METHOD clc; clearall; x(1)=input('Enter the value x0 = '); xn=input('Enter the value of xn = '); y(1)=input('Enter the value of y0 = '); f=input('Enter the function of yd, f = '); h=input('Enter the value of h = '); ac=input('Enter the value of accuracy, ac = '); n=((xn-x(1))/h); for i=1:n+1 x(i+1)=x(i)+h; y(i+1)=y(i)+h*f(x(i),y(i)); e=1; while (abs (e)>ac) yc(i+1)=y(i)+(h/2)*(f(x(i),y(i))+(f(x(i)+h,y(i+1)))); e=yc(i+1)-y(i+1); y(i+1)=yc(i+1); end ; fprintf('\n At x=%f y=%f',x(i+1),y(i+1)); end ;
Output: Enter the value of x0 = 1 Enter the value of xn = 1.2 Enter the value of y0 = 2.2 Enter the function of yd, f = @(x,y)(sqrt(x+y)) Enter the value of h = 0.1 Enter the value of accuracy, ac = 0.0001 At x=1.100000 y=2.382753 At x=1.200000 y=2.573186
Solver: Save .m file with name MODEULRsolver functiondy=f(x,y) dy= sqrt(x+y);
y =
2.2000 2.3828 2.5732
% PROGRAM FOR RUNGE KUTTA 2ND ORDER clc; clear all; x(1)=input('Enter the value of x0 = '); xn=input('Enter the value of xn = '); y(1)=input('Enter the value of y0 = '); f=input('Enter the function of yd , f = '); h=input('Enter the value of h = '); n=((xn-x(1))/h); for i=1:n+1 x(i+1)=x(i)+h; k1=h*f(x(i),y(i)); k2=h*f(x(i)+h,y(i)+k1); dy=(k1+k2)/2; y(i+1)=y(i)+dy; fprintf('\n At x = %f,
y = %f',x(i+1),y(i+1));
end ;
Output: Enter the value of x0 = 0 Enter the value of xn = 0.4 Enter the value of y0 = 1.2 Enter the function of yd , f = @(x,y)((x^2)/(2*y)) Enter the value of h = 0.2 At x = 0.200000, y = 1.201667 At x = 0.400000, y = 1.209970
Solver: Save .m file with name RK2solver
functiondy=f(x,y) dy=(x^2)/(2*y);
@ command window: >> [x,y]=ode45('RK2solver',[0:0.2:0.4],1.2)
% Program for parabolic explicit method clc; clearall; x(1)=input('Enter xo= '); xn=input('Enter xn= '); t(1)=input('Enter to= '); tm=input('Enter tm= '); c=input('Enter c Square= '); h=input('Enter increment in x,h= '); k=input('Enter increment in t,k= '); p=input('Enter value for u(0,t) i.e first column= '); q=input('Enter value for u(n,t) i.e last column= '); r=input('Enter function for u(x,0) for 0
a=k*c/(h^2); n=(xn-x(1))/h; m=(tm-t(1))/k;
for i=1:n x(i+1)=x(i)+h; end ;
for i=1:m t(i+1)=t(i)+k; end ;
for j=1:m+1 u(1,j)=p; u(n+1,j)=q; end ; for i=2:n u(i,1)=r(x(i),t(1)); end ;
for j=1:m
Output: Enter xo= 0 Enter xn= 0.5 Enter to= 0 Enter tm= 0.03 Enter c Square= 1 Enter increment in x,h= 0.1 Enter increment in t,k= 0.01 Enter value for u(0,t) i.e first column= 1 Enter value for u(n,t) i.e last column= 1 Enter function for u(x,0) for 0
the table for u is as under: 1.0000
1.0000
1.0000
1.0000
1.2000
1.2000
1.2000
1.2000
1.4000
1.4000
1.4000
0.4000
1.6000
1.6000
0.6000
2.6000
1.8000
0.8000
1.8000
-0.2000
1.0000
1.0000
1.0000
1.0000
% Program for Hyperbolic equation clc; clearall;
cs=input('Enter c Square= '); x(1)=input('Enter xo= '); xn=input('Enter xn= '); t(1)=input('Enter to= '); h=input('Enter increment in x,h= '); k=input('Enter increment in t,k= '); p=input('Enter value for u(0,t) i.e first column= '); q=input('Enter value for u(xn,t) i.e last column= '); r=input('Enter function for u(x,0) for 0
