Runge-Kutta Method

%% Write a MATLAB program to find y for x =0.2 given that dy/dx =

%% x+y and y=1 when x=0 using Runga-Kutta fourth order method.



% Given

x0 = 0;

y0 = 1;

%let

h = 0.2;

fxy = x0+y0;

k1 =h*fxy;

x0=x0+h/2;y0=y0+k1/2;

fxy = x0+y0;

k2=h*fxy;

y0=y0-k1/2+k2/2;

fxy = x0+y0;

k3=h*fxy;

x0=x0+h/2;y0=y0-k2/2+k3;

fxy = x0+y0;

k4=h*fxy;

k=(k1+2*k2+2*k3+k4)/6

%then y = y0+k

fprintf (' Approximate value of y for x = %2.3f is = %2.4f \n',x0,1+k);

Euler's Method

%% Write a MATLAB program to find y for x =0.1 given that dy/dx =

%% (y-x)/(y+x)taking h=0.02 using Euler's Method.



% Given

x0 = 0;

y0 = 1;

h = 0.02;

x =0.1;

x1 =x0;

y1=y0;

i=0;

fprintf (' x%1.0f = %2.2f, y%1.0f = %2.4f \n',i,x1,i,y1);

while x1<x

y1= y1+h*(y1-x1)/(y1+x1);

x1=x1+h;

i=i+1;

fprintf (' x%1.0f = %2.2f, y%1.0f = %2.4f \n',i,x1,i,y1);

end

Integration by Gaussian Integration Method

%% Write a MATLAB program to integrate the equation f(x) = dx/(1+x) within

%% the limits of 0 to 1 with three points using Gaussian Integration Method.

% Given a= 0;

b = 1;

n = 3;
%% writing x in terms of u ( x = a1*u + b1)with limits of -1 to 1a1 = (b-a)/2;

b1 = (b+a)/2;

dx = a1;

%a1=1; b1=0;

% Using the table for w and x values for n=3

u1 = -0.77460; u2 = 0.0; u3 = 0.77460;

w1 = 0.55555; w2 = 0.88889; w3 = 0.55555;

% fu1 = du/1+(u1)^2 where dx=a1*du

fu1 = a1/(1+(a1*u1+b1));

fu2 = a1/(1+(a1*u2+b1));

fu3 = a1/(1+(a1*u3+ b1));



I = w1*fu1+w2*fu2+w3*fu3;

fprintf (' The integral by Gaussian integration method is = %2.5f \n',I);

Integration by Simpsons 1/3 rule

%% Write a MATLAB program to integrate the equation f(x) = dx/(1+x^2) within

%% the limits of 0 to 6 by taking 6 subintervals by using simpsons 1/3 rule.



% Given

x0 = 0;

xn = 6;

n = 6;

h = (xn-x0)/n;

y0 = 1/(1+x0^2);

yn = 1/(1+xn^2);

% sum of integral terms is given by y0+yn+4*(y1+y3+....yn-1)+2*(y2+y4+...Yn-2)

s = 0;

for i = 1:2:n-1

s = s + 4*1/(1+(x0+i*h)^2);

end

for i = 2:2:n-1

s = s + 2*1/(1+(x0+i*h)^2);

end

s = s+y0+yn;

% integral is given by h/3(sum of integral terms)

integral = (h/3)*s;

fprintf (' The integral by simpsons 1/3 rule is = %2.5f \n',integral);

MATLAB code to find inverse of a Matrix by Gauss Emilination

  %% Write a program which can calculate inverse of a given matrix by

  %% by Gauss Elimination Method using complete pivot.  

    M = [1 2 13;3 24 5;25 65 7]; % Given matrix whose inverse is to be found

    A = M;

    Ni = length(A(1,:)); % number of rows of given matrix

    Nj =length(A(:,1));   

    if abs(det(A))<0.00000000001

        fprintf('The matrix determinant is zero, the result found may not be correct');

    end     

    px =1:Ni;

    qy =1:Ni;  
 

    %%main loop; find pivot entry in columns qy[k], qy[k+1], ...,qy[n-1]   

    %% and in rows px[k], px[k+1],.....px[n-1]

  

    for k = 1: Ni-1

        pct = k; qdt = k;

        aet =0;

        for i=k:Ni

            for j=k:Ni

                tmp = abs(A(uint8(px(i)),uint8 (qy(j))));

                if (tmp>aet)

                    aet = tmp;  pct = i;qdt = j;

                end

            end

        end

      

        if ~aet

            fprintf('pivot is zero, No inverse possible \n');

        end

     

     % Swap the rows and columns to get highest number in the diagonal

        px([k pct]) = px([pct k]);%//swap px[k] and px[pct]                        

        qy([k qdt]) = qy([qdt k]) ; 

        %% eliminate column entries (make them zero) logically below and right to the entry mx[px[k]][qy[k]]

        for i=k+1:Ni   

            if A(px(i),qy(k)) ~= 0  % if the entry is not equal to zero        

                mult = A(px(i),qy(k))/A(px(k),qy(k)); % get the multiplying factor

                A(px(i),qy(k))= mult;

                for j=k+1:Ni

                    A(px(i),qy(j)) = A(px(i),qy(j)) - mult*A(px(k),qy(j)); % entry will become zero

                end

            end

        end

    end

  

    AI = eye(size(A));  %% Identiy Matrix

    %%forward substitution 

    for k=1:Ni

        for i =2:Ni

            for j =1:i-1

                AI(px(i),k) = AI(px(i),k) - A(px(i),qy(j))*AI(px(j),k);

            end

        end

    end      

    InMat = zeros(Ni,Nj); %InMat will be the required inverse after back substitution.  

    %% backward substitution

    for k=1:Ni

        for i = Ni:-1:1

            for j = i+1:Ni

                AI(px(i),k) = AI(px(i),k) - A(px(i),qy(j))* InMat(qy(j),k);

            end

        InMat(qy(i),k) = AI(px(i),k) / A(px(i),qy(i));

        end

    end

    %% wee can also find error of our result by subtracting the identity matrix

    %% from the product of given matrix (M) and our result (InMat)

    I = eye(Ni,Nj);

    error = InMat*M - I;

   %% output of the results

    fprintf('The given matrix is \n\n');

    fprintf(' %f %f %f \n',M);

    fprintf('The inverse of the matrix is \n\n');

    fprintf(' %f  %f  %f \n',InMat);

    fprintf('The error in the calculate of the invese of matrix is \n\n');

    fprintf('\n %f  %f  %f \n',error);