Matlab Codes for Numerical Methods

Here are the Matlab codes for solving few Numerical Methods problems. As per our course requirement for Mechanical Engineering, sophomore course, we had to study Numerical Methods and the Matlab programming was in practical portion. And I have tried to share few codes with you for your convenience with codes.

This page includes Bisection Method, Regula-Falsi Method, Secant Method, General Newton Method, Rapshon Newton Method, Euler Method, Jacobi Method, Gauss Seidal Method, Simpson(1/3), Simpson(3/8) and Trapeziodal Rule. 
I have tried all the codes by making function file.

To execute the program first of all you have to define all the required function and variables in command window. For example to execute the program codes for Bisection Method, do the following in command window :-

>> f =inline('x^2+3*x+10'); (put your own function )

>>a=0.1; (put your own value)

>>b=2; (put your own value )

>>tol=0.0001; (put your own value )

>> bisection(f, a, b, tol) and press enter

OR you can directly use

>>bisection(f,0.1,2,0.000.1)

The main thing to remember is to satisfy the required function prototype.

Bisection Method

function[]=bisection(f,a,b,tol);

n=round(log(abs(b-a)/tol)/log(2));    %Calcuation_of_no_of_iterations

i=1;

if f(a)*f(b)>0

   disp('error');

   end

   fprintf('iteration\ta\t\tb\t\tp\t\tf(p)\n');     %Display_table

while i<n

   p=(a+b)/2;

   fprintf('%d\t\t%.2f\t\t%.2f\t\t%.2f\t\t%.2f\n\n\n',i,a,b,p,f(p));      %Display_table _with_successive _iterations

   if f(a)*f(p)<0    &Condition_check

      a=a;

      b=p;

   else

      a=b;

      b=p;

   end;

   i=i+1;

end;

   fprintf('The solution is %.4f\n', p);

Regula-Falsi Method

function[]=falsi(f,a,b,tol);

i=1;

if f(a)*f(b)>0

    disp('error');

end;

while i<50

   

    x=(a*f(b)-b*f(a))/(f(b)-f(a));

    if f(a)*f(x)>0

        a=a;

        b=x;

    else

        a=x;

        a=b;

        if abs(x-a)||abs(x-b)==tol

        fprintf('\nthe solution is %0.4f\n',x);

        end;

        break;

    end;

        i=i+1;

    end;

      

Secant Method

function[]=secant(f,a,b,tol);

i=1;

if f(a)*f(b)>0

    disp('error');

end;

while i<50

   

    x=(a*f(b)-b*f(a))/(f(b)-f(a));

  

   

        a=b;

        b=x;

        if abs(x-a)||abs(x-b)==tol

        fprintf('\nthe solution is %0.4f\n',x);

      

        break;

    end;

        i=i+1;

    end;

General Newton Method

function[]=newtongeneral(f,g,x0,p,tol);

x(1)=x0;

i=1;

while i<50         %Taking_arbitrary_no_50__but_You_can_also_use_infinite_loop_using_while(1)

    x(i+1)=x(i)-p*(f(x(i))/g(x(i)));

    if abs(x(i+1)-x(i))<=tol

        X=x(i+1);

        fprintf('The solution is x = %.4f\n',X);

        break;

    else

        i=i+1;

    end;

end;

       

Raphson Newton Method

function[]=newtonrap(f,g,x0,tol);

x(1)=x0;

i=1;

while (1)

    x(i+1)=x(i)-(f(x(i)))/(g(x(i)));

    if abs(x(i+1)-x(i))<=tol

        X=x(i+1);

        fprintf('The solution is x=%.4f\n',X);

        break;

    else

       i=i+1;

    end;

end;

   

Euler Method

function[]=euler(f,x0,y0,h);

i=1;

x(1)=x0;

y(1)=y0;

n=(y0-x0)/h;

while i<n

    x(i)=x0+i*h;

    y(i+1)=y(i)+h*f(x(i),y(i));

    i=i+1;

end;

fprintf('The solutions are %.4f %.4f %.4f', y(1),y(2),y(3));  %_You_can_use_y(n)_upto_n_values_

Jacobi Method

function[]=jacobi(f,g,h,e);

i=1;

A(1)=f(0,0);

B(1)=g(0,0);

C(1)=h(0,0);

while i<50

    A(i+1)=f(B(i),C(i));

    B(i+1)=g(C(i),A(i));

    C(i+1)=h(A(i),B(i));

   

    if (abs(A(i+1)-A(i))==e && abs(B(i-1)-B(i))==e && abs(C(i-1)-C(i))==e)

        X=A(i+1);

        Y=B(i+1);

        Z=c(i+1);

        break;

    else

     i=i+1;    

    end;

    end;

fprintf('\nThe solutions are X= %0.4f\n Y= %0.4f \n Z= %0.4f\n',A(i),B(i),C(i));

Gauss Seidal Method

function[]=gaussseidal(f,g,h,e);

i=1;

A(1)=f(0,0);

B(1)=g(0,0);

C(1)=h(0,0);

while (1)

    A(i+1)=f(B(i),C(i));

    B(i+1)=g(C(i),A(i+1));

    C(i+1)=h(A(i+1),B(i+1));

   

    if (abs(A(i+1)-A(i))==e && abs(B(i-1)-B(i))==e && abs(C(i-1)-C(i))==e)

        break;

        X=A(i+1);

        Y=B(i+1);

        Z=c(i+1);

       

    else

     i=i+1;    

    end;

    end;

fprintf('\nThe solutions are X= %0.4f\n Y= %0.4f \n Z= %0.4f\n',A(i),B(i),C(i));

Simpson Rule (1/3)

function[]=simpsonthird(f,a,b,h);

n=(b-a)/h;

y0=f(a);

yn=f(b);

i=1;

sum1=0;

while (2*i-1)<n

    x(2*i-1)=a+(2*i-1)*h;

    y(2*i-1)=f(x(2*i-1));

    sum1=sum1+y(2*i-1);

    i=i+1;

end;

i=1;

sum2=0;

while 2*i<n-1

    x(2*i)=a+(2*i)*h;

    y(2*i)=f(x(2*i));

    sum2=sum2+y(2*i);

    i=i+1;

end;

I=(h/3)*(y0+yn+4*sum1+2*sum2);

fprintf('\nThe integrated value is %0.4f\n',I);

Simpson Rule (3/8)

function[]=simpsoneight(f,a,b,h);

n=(b-a)/h;

y0=f(a);

yn=f(b);

i=1;

j=1;

k=1;

sum1=0;

sum2=0;

sum3=0;

while (3*i)<n-2

    x(3*i)=a+(3*i)*h;

    y(3*i)=f(x(3*i));

    sum1=sum1+y(3*i);

    i=i+1;

end;

while (3*j-2)<n

    x(3*j-2)=a+(3*j-2)*h;

    y(3*j-2)=f(x(3*j-2));

    sum2=sum2+y(3*j-2);

    j=j+1;

end;

while (3*k-1)<n-1

    x(3*k-1)=a+(3*k-1)*h;

    y(3*k-1)=f(x(3*k-1));

    sum3=sum3+y(3*k-1);

    k=k+1;

end;

I=(3*h/8)*(y0+yn+2*(sum1)+3*(sum2+sum3));

fprintf('\nThe Integrated value is %0.4f\n',I);

Trapezoidal Rule

function []=trapezoidal(f,a,b,h);

n=(b-a)/h;

y0=f(a);

yn=f(b);

i=1;

sum=0;

while i<n

    x(i)=a+i*h;

    y(i)    =f(x(i));

    sum=sum+y(i);

    i=i+1;

end;

I=(h/2)*(y0+yn+2*sum);

fprintf('\nThe Integrated value is %0.4f\n\n',I);

Thank you for looking this page. Please comment. Have a good day.

Read more »