MATLAB code for Conjugate gradient method

function [x]= conjugate_gradient_method(A,b,x)

    r0=b-A*x;% A is a square matrix of nxn and x is a vector of nx1

             % b is also a vector of nx1 order

             % equation is of form AX=B

    d=r0;

    for i= 1:length(b)

        alpha = (r0'*r0)/(d'*A*d);fprintf('alpha = %f\n\n',alpha);

        x = x + alpha * d;fprintf('x = %f\n',x);fprintf('\n');

        r1 = r0 - alpha * A * d;fprintf('r = %f\n',r1);fprintf('\n');

        if sqrt(r1'*r1)< 1e-10

            break;

        end

        beta = (r1'*r1)/(r0'*r0);fprintf('beta = %f\n',beta);fprintf('\n');

        d = r1 + beta * d;fprintf('d = %f\n',d);fprintf('\n');

        r0 = r1;

        fprintf('----------------------------------\n\n');

    end

end

Read More

MATLAB Programming code for Gauss-siedal method

n=input('Enter number of simultaneous equation : ');

a=input('Enter the coefficient matrix [a] : ');

b=input('Enter the constant column matrix {b} : ');

x=input('Initialze all values {x} : ');

for c=1:100

for i = 1:n

    sum=0;

    for k = 1:n

        if k~=i

            sum=sum+a(i,k)*x(k);

        end

    end

    x(i)=(b(i)-sum)/a(i,i); %suppose i=1 then x(1) gives value of 1st row

end

end

disp(x) %disp display each value in new line

   

Read More

Code for Gauss-siedal method in C programming

#include<stdio.h>
void main()
{
    float a[50][50],x[50],b[50],sum,e[50],y[50],err;
    int n,i,j,k,itr,c,r=0;
    printf("Enter order of matrix, no. of iteration and error : ");
    scanf("%d%d%f",&n,&itr,&err);
    for(i=1;i<=n;i++)
    {
        for(j=1;j<=n;j++)
        {
            printf("Enter a[%d][%d] = ",i,j);
            scanf("%f",&a[i][j]);
        }
    }
    for (i=1;i<=n;i++)
    {
        printf("Enter b[%d] = ",i);
        scanf("%f",&b[i]);
    }
    for (i=1;i<=n;i++)
    {
        printf("Enter x[%d] = ",i);
        scanf("%f",&x[i]);
        y[i]=x[i];
    }


        for(c=1;c<=itr;c++)
        {
            printf("\niteration %d",c);
            for(i=1;i<=n;i++)
            {

            sum=0;
            for(k=1;k<=n;k++)
            {
                if(k==i)
                    continue;
                    sum=sum+a[i][k]*x[k];

            }
            x[i]=(b[i]-sum)/a[i][i];
            printf("\tx%d=%f\t",i,x[i]);
            e[i]=fabs(y[i]-x[i]);
            //printf("\te%d=%f\t",i,e[i]);
            y[i]=x[i];
            }
            for(j=1;j<=n;j++)
            {
                if(e[j]<=err)
                {
                    r++;
                }
                if(r==n)
                break;
            }
            if(r==n)
                break;

        }
    if(r==n)
    {
        for(i=1;i<=n;i++)
        {
            printf("\n\nx%d=%f",i,x[i]);
        }
    }
    else
    {
        printf("\nSolution couldn't be converged !");
    }

         getch();
}

Read More

Program code in MATLAB to calculate the displacement of axially loaded tapered bar fixed at one end by Finite Element Method

                                                             

                                                              

Program code in MATLAB 

function  tappered_bar(At,Ab,l,fb,E,n,filename)

    y0=At;

    K=zeros(n);k=linspace(0,0,n);f=linspace(0,0,n)';f(n)=fb; )%At=area at fixed end Ab=area at free end

    op=fopen(filename,'wt'); %fb=axial force at free end and n=no. of element

fprintf(op,'=======================================================================\n');

fprintf('=======================================================================\n');

fprintf(op,'\t\t\tTappered bar solution by FEM\n'); fprintf('\t\t\t\t\tTappered bar solution by FEM\n');

fprintf(op,'-----------------------------------------------------------------------\n');

fprintf('-----------------------------------------------------------------------\n');

fprintf(op,'Maximum area = %f\t\t\t',At); fprintf(op,'Minimum area = %f\n',Ab);

fprintf('Maximum area = %f\t\t\t\t',At); fprintf('Minimum area = %f\n',Ab);

fprintf(op,'Force at minimum area = %f\t',fb); fprintf(op,'Lenght = %f\n',l);

fprintf('Force at minimum area = %f\t',fb); fprintf('Lenght = %f\n',l);

fprintf(op,'Modulus of elasticity = %f\t',E); fprintf(op,'No. of element = %f\n',n);

fprintf('Modulus of elasticity = %f\t',E); fprintf('No. of element = %f\n',n);

fprintf(op,'-----------------------------------------------------------------------\n');

fprintf('-----------------------------------------------------------------------\n');

    for i = 1:n

        yi = (Ab-At)*i/n+At;

        Ai= (y0 +yi)/2;

        k(i)=n*E*Ai/l;

        y0=yi;

        fprintf(op,'A%d = %f\t',i,Ai);fprintf('A%d = %f\t',i,Ai);

        fprintf(op,'k%d = %f\n',i,k(i));fprintf('k%d = %f\n',i,k(i));

    end

fprintf(op,'-----------------------------------------------------------------------\n');

fprintf('-----------------------------------------------------------------------\n');

    for i=1:n

        for j=1:n

            if abs(i-j)==1

                if i>j

                    K(i,j)=-k(i);

                elseif j>i

                        K(i,j)=-k(j);

                end

            elseif i==n && j==n

                    K(i,j)=k(i);

            elseif i~=n && j~=n && i==j

                    K(i,j)=k(i)+k(i+1);

            end

        end   

    end

    u =  inv(K)*f;

    disp('K = ');fprintf(op,'K = \n');

    disp(K);

    for i=1:n

        for j=1:n

            fprintf(op,'%f\t',K(i,j));

        end

        fprintf(op,'\n');

    end

fprintf(op,'-----------------------------------------------------------------------\n');

fprintf('-----------------------------------------------------------------------\n');

    for i=1:n

        fprintf(op,'u%d = %f\n',i,u(i));fprintf('u%d = %f\n',i,u(i));

        %u_approx=ui;

    end

fprintf(op,'=======================================================================\n');

fprintf('=======================================================================\n');

fprintf(op,'\t\t\tExact solution of Tappered bar\n');fprintf('\t\t\t\t\tExact solution of Tappered bar\n');

fprintf(op,'-----------------------------------------------------------------------\n');

fprintf('-----------------------------------------------------------------------\n');

    syms x;

    Area_at_x = (Ab-At)*x/l+At;

    del= (fb/E)*int(1/Area_at_x,x,0,l);

    fprintf(op,'\t\t\t\tdel = %f\n',del);fprintf('\t\t\t\t\t\t\tdel = %f\n',del);

fprintf(op,'=======================================================================\n');

fprintf('=======================================================================\n');

    plot(n,u(n),'r*',n,del,'gx');

    xlabel('no. of Elements for FEM solution');ylabel('Displacement');title('Comapring FEM and Exact solution');

    legend('FEM solution','Exact solution');

    fclose(op);

end

OUTPUT when n= 3

=======================================================================

                                                Tappered bar solution by FEM

--------------------------------------------------------------------------------------------------------------------

Maximum area = 0.031416                        Minimum area = 0.007854

Force at minimum area = 10.000000           Lenght = 2.000000

Modulus of elasticity = 1.000000                     No. of element = 3.000000

-------------------------------------------------------------------------------------------------------------------

A1 = 0.027489          k1 = 0.041233

A2 = 0.019635          k2 = 0.029452

A3 = 0.011781          k3 = 0.017671

-------------------------------------------------------------------------------------------------------------------

K =

0.070686  -0.029452 0.000000 

-0.029452 0.047124  -0.017671

0.000000  -0.017671 0.017671 

-------------------------------------------------------------------------------------------------------------------

u1 = 242.521818

u2 = 582.052363

u3 = 1147.936605

=======================================================================

                                                Exact solution of Tappered bar

-------------------------------------------------------------------------------------------------------------------

                                                                del = 1176.723201

=======================================================================

OUTPUT when n=4

=======================================================================

                                                Tappered bar solution by FEM

--------------------------------------------------------------------------------------------------------------------

Maximum area = 0.031416                        Minimum area = 0.007854

Force at minimum area = 10.000000           Lenght = 2.000000

Modulus of elasticity = 1.000000                  No. of element = 4.000000

--------------------------------------------------------------------------------------------------------------------

A1 = 0.028471          k1 = 0.056941

A2 = 0.022580          k2 = 0.045160

A3 = 0.016690          k3 = 0.033379

A4 = 0.010799          k4 = 0.021598

--------------------------------------------------------------------------------------------------------------------

K =

0.102102  -0.045160 0.000000  0.000000 

-0.045160 0.078540  -0.033379 0.000000 

0.000000  -0.033379 0.054978  -0.021598

0.000000  0.000000  -0.021598 0.021598 

--------------------------------------------------------------------------------------------------------------------

u1 = 175.619248

u2 = 397.052212

u3 = 696.637987

u4 = 1159.634185

=======================================================================

                                                Exact solution of Tappered bar

-------------------------------------------------------------------------------------------------------------------

                                                                del = 1176.723201

=======================================================================

Read More