Doolittle Algorithm: LU Decomposition Algorithm

Arogya Thapa Magar at March 29, 2019 0

DOOLITTLE Algorithm

In the numerical method Doolittle Algorithm : LU Decomposition Algorithm (where LU stands for Lower and Upper and also called LU factorization Algorithm) factors a matrix as the product of a lower triangular matrix and an upper triangular matrix. Computers usually solve square systems of linear equations using the LU decomposition, and it is also a key step when inverting a matrix, or computing the determinant of a matrix.

Let A be a square matrix. An LU factorization algorithm refers to the factorization of A, with proper row and/or column orderings or permutations, into two factors, a lower triangular matrix L and an upper triangular matrix U, A=LU.
Assume that A has a Crout factorization A = LU.

Doolittle Algorithm :-

It is always possible to factor a square matrix into a lower triangular matrix and an upper triangular matrix. That is, [A] = [L][U]

Doolittle’s method provides an alternative way to factor A into an LU decomposition without going through the hassle of Gaussian Elimination.

For a general n×n matrix A, we assume that an LU decomposition exists, and write the form of L and U explicitly. We then systematically solve for the entries in L and U from the equations that result from the multiplications necessary for A=LU.

#include<stdio.h>
#include<conio.h>
#include<math.h>

int main(){
 int n,i,k,j;
 float sum=0,a[10][10],u[10][10],l[10][10],b[10],x[10],z[10];
 
 printf("Do-Little LU Decomposition");
 printf("\nEnter Dimension of System of equation\n");
 scanf("%d",&n);
 printf("\nEnter the coefficients of the Matrix\n");
 for(i=0;i<n; i++)
 for(j=0;j<n;j++){
  scanf("%f",&a[i][j]);
 }
 printf("Enter RHS vector\n");
 for(i=0; i<n; i++){
  scanf("%f",&b[i]);
 }
 
 //Compute Elements of L and U matrix
 for(j=0; j<n; j++)
  u[0][j]=a[0][j];
 
 for(i=0;i<n;i++)
  l[i][i]=1;
 
 for(i=1; i<n; i++)
  l[i][0]=a[i][0]/u[0][0];
 
 
 for(j=1;j<n;j++){
  for(i=1; i<=j;i++){
   for(k=0;k<=i-1;k++){
    sum=sum+(l[i][k]*u[k][j]);
   }
   u[i][j]=a[i][j]-sum;
   sum=0;
  }
  
  for(i=j+1;i<n;i++){
   for(k=0;k<=j-1;k++){
    sum=sum+(l[i][k]*u[k][j]);
   }
   l[i][j]=(a[i][j]-sum)/u[j][j];
   sum=0;
  }
 }
 
 //Solve for Z using forward subsitution
 z[0]=b[0];
 for(i=1;i<n;i++){
  for(j=0; j<i; j++)
  sum=sum+(l[i][j]*z[j]);
  z[i]=b[i]-sum;
  sum=0;
 }
 // solve for X using Backward subsitution
 x[n-1]= z[n-1]/u[n-1][n-1];
 for(i=n-2;i>=0;i--){
  for(j=i+1;j<n;j++)
  sum=sum+(u[i][j]*x[j]);
  x[i]=(z[i]-sum)/u[i][i];
  sum=0;
 }
 printf("\nSolution:\n");
 for(i=0; i<n;i++){
  printf("x%d=%f\t",i+1,x[i]);
 }
 getch();
 return 0;
}

Numerical Methods

Matrix Inverse : Matrix Inversion Technique

Arogya Thapa Magar at March 26, 2019 0

Matrix Inverse : Matrix Inversion Technique

The inverse of a matrix

The inverse of a square n × n matrix A is another n × n matrix denoted by A^-1 such that

AA^-1=A^-1A=I

where I is the n × n identity matrix. That is, multiplying a matrix by its inverse produces an identity matrix. Not all square matrices have an inverse matrix. If the determinant of the matrix is zero, then it will not have an inverse, and the matrix is said to be singular. Only non-singular matrices have inverses.

A formula for finding the inverse

Given any non-singular matrix A, its inverse can be found from the formula A^-1 = adj A |A| where adj A is the adjoint

matrix and |A| is the determinant of A. The procedure for finding the adjoint matrix is given below.

Finding the adjoint matrix

The adjoint of a matrix A is found in stages:

Find the transpose of A, which is denoted by A^T. The transpose is found by interchanging the rows and columns of A. So, for example, the first column of A is the first row of the transposed matrix; the second column of A is the second row of the transposed matrix, and so on.

The minor of any element is found by covering up the elements in its row and column and finding the determinant of the remaining matrix. By replacing each element of A^T by its minor, we can write down a matrix of minors of A^T.

The cofactor of any element is found by taking its minor and imposing a place sign according to the following rule

This means, for example, that to find the cofactor of an element in the first row, second column, the sign of the minor is changed. On the other hand to find the cofactor of an element in the second row, second column, the sign of the minor is unaltered. This is equivalent to multiplying the minor by ‘+1’ or ‘−1’ depending upon its position. In this way, we can form a matrix of cofactors of A^T. This matrix is called the adjoint of A, denoted adj A. The matrix of cofactors of the transpose of A, is called the adjoint Matrix, adj A.

C program to find value using matrix inversion technique

#include<stdio.h>
#include<conio.h>
#include<math.h>

int main(){
 int n,i,k,j,p,q;
 float pivot,term,a[10][10];
 printf("Matrix Inversion");
 printf("\nEnter Dimension of System of equation\n");
 scanf("%d",&n);
 printf("\nEnter the coefficients of the Matrix\n");
 for(i=0;i<n; i++)
 for(j=0;j<n;j++){
  scanf("%f",&a[i][j]);
 }
 
 for(i=0; i<n; i++){
  for(j=n;j<2*n;j++){
   if(i==j-n)
    a[i][j]=1;
   else
    a[i][j]=0;
  }
 }
 
 for(k=0; k<n; k++){
  pivot=a[k][k];
  for(p=0; p<2*n;p++)
  a[k][p]=a[k][p]/pivot;
  for(i=0;i<n;i++){
   term=a[i][k];
   if(k!=i)
   for(j=0;j<2*n;j++){
    a[i][j]=a[i][j]-a[k][j]*term;
   }
  }
 }
 
 printf("\nMatrix Inverse is:\n");
 for(i=0; i<n;i++){
  for(j=n;j<2*n;j++)
  printf("%f\t",a[i][j]);
  printf("\n");
 }
 getch();
return 0;
}

Java code to find inverse of matrix

//This is sample program to find the inverse of a matrix
 
import java.util.Scanner;
 
public class Inverse 
{
 public static void main(String argv[]) 
 {
 Scanner input = new Scanner(System.in);
 System.out.println("Enter the dimension of square matrix: ");
 int n = input.nextInt();
 double a[][]= new double[n][n];
 System.out.println("Enter the elements of matrix: ");
 for(int i=0; i<n; i++)
 for(int j=0; j<n; j++)
 a[i][j] = input.nextDouble();
 
 double d[][] = invert(a);
 
 System.out.println("The inverse is: ");
 for (int i=0; i<n; ++i) 
 {
 for (int j=0; j<n; ++j)
 {
 System.out.print(d[i][j]+" ");
 }
 System.out.println();
 }
 input.close();
 } 
 
 public static double[][] invert(double a[][]) 
 {
 int n = a.length;
 double x[][] = new double[n][n];
 double b[][] = new double[n][n];
 int index[] = new int[n];
 for (int i=0; i<n; ++i) 
 b[i][i] = 1;
 
 // Transform the matrix into an upper triangle
 gaussian(a, index);
 
 // Update the matrix b[i][j] with the ratios stored
 for (int i=0; i<n-1; ++i)
 for (int j=i+1; j<n; ++j)
 for (int k=0; k<n; ++k)
 b[index[j]][k]
 -= a[index[j]][i]*b[index[i]][k];
 
 // Perform backward substitutions
 for (int i=0; i<n; ++i) 
 {
 x[n-1][i] = b[index[n-1]][i]/a[index[n-1]][n-1];
 for (int j=n-2; j>=0; --j) 
 {
 x[j][i] = b[index[j]][i];
 for (int k=j+1; k<n; ++k) 
 {
 x[j][i] -= a[index[j]][k]*x[k][i];
 }
 x[j][i] /= a[index[j]][j];
 }
 }
 return x;
 }
 
// Method to carry out the partial-pivoting Gaussian
// elimination. Here index[] stores pivoting order.
 
 public static void gaussian(double a[][], int index[]) 
 {
 int n = index.length;
 double c[] = new double[n];
 
 // Initialize the index
 for (int i=0; i<n; ++i) 
 index[i] = i;
 
 // Find the rescaling factors, one from each row
 for (int i=0; i<n; ++i) 
 {
 double c1 = 0;
 for (int j=0; j<n; ++j) 
 {
 double c0 = Math.abs(a[i][j]);
 if (c0 > c1) c1 = c0;
 }
 c[i] = c1;
 }
 
 // Search the pivoting element from each column
 int k = 0;
 for (int j=0; j<n-1; ++j) 
 {
 double pi1 = 0;
 for (int i=j; i<n; ++i) 
 {
 double pi0 = Math.abs(a[index[i]][j]);
 pi0 /= c[index[i]];
 if (pi0 > pi1) 
 {
 pi1 = pi0;
 k = i;
 }
 }
 
 // Interchange rows according to the pivoting order
 int itmp = index[j];
 index[j] = index[k];
 index[k] = itmp;
 for (int i=j+1; i<n; ++i) 
 {
 double pj = a[index[i]][j]/a[index[j]][j];
 
 // Record pivoting ratios below the diagonal
 a[index[i]][j] = pj;
 
 // Modify other elements accordingly
 for (int l=j+1; l<n; ++l)
 a[index[i]][l] -= pj*a[index[j]][l];
 }
 }
 }
}

Numerical Methods

Improved Euler's Method and Euler's Method with their C Program

Arogya Thapa Magar at March 24, 2019 0

Improved Euler’s Method

The Improved Euler’s method, also known as the Heun formula or the average slope method, gives a more accurate approximation than the Euler rule and gives an explicit formula for computing y_n+1. The basic idea is to correct some errors of the original Euler's method. The syntax of the Improved Euler’s method is similar to that of the trapezoid rule, but the y value of the function in terms of y_n+1 consists of the sum of the y value and the product of h and the function in terms of x_n and y_n.

Improved Euler formula or the average slope method is commonly referred to as the Heun method:

yn+1=yn

+h/2[f(xn,yn)+f(xn+1,yn+hf(xn,yn))],n=0,1,2,….

Since it is actually the simplest version of the predictor-corrector method, the recurrence can be written as

pn+1

yn+1=yn+hf(xn,yn),

yn+1

=yn+h/2[f(xn,yn)+f(xn+1,pn+1)],n=0,1,2,….

C program to find the value using Heun's method

#include<stdio.h>
#include<math.h>
#include<conio.h>
#define f(x,y) 2*(x)+(y)

int main(){
 float x,xp,x0,y0,y,h,m1,m2;
 printf("\n Heuns Method\n");
 printf("\nEnter initial values of x and y\n");
 scanf("%f%f",&x0,&y0);
 printf("\nEnter x at which function to be evaluated\n");
 scanf("%f",&xp);
 printf("\nEnter the step size \n");
 scanf("%f",&h);
 y=y0;
 x=x0;
 for(x=x0;x<xp;x=x+h){
  m1=f(x,y);
  printf("\nm1=%f",m1);
  m2=f(x+h,y+(h*m1));
  printf("\nm2=%f",m2);
  y=y+h/2*(m1+m2);
  printf("\ny=%f",y);
 }
 
 printf("Function value at x=%f is %f\n",xp,y);
 getch();
 return 0; 
 
}

Euler's Method

Euler's Method assumes our solution is written in the form of a Taylor’s Series.

That is, we'll have a function of the form:

This gives us a reasonably good approximation if we take plenty of terms, and if the value of

For Euler's Method, we just take the first 2 terms only.

The last term is just

h times our dx/dy

expression, so we can write Euler's Method as follows:

How do we use this formula?

We start with some known value for y, which we could call y0. It has this value when X

=x

. (We make use of the initial value

(

x0,Y0)).

The result of using this formula is the value for y, one h step to the right of the current value. Let's call it

Y

1So we have:

where
$Y$ 1 is the next estimated solution value;
Y0 is the current value;
h is the interval between steps; and
$f$ (x0,y0) is the value of the derivative at the starting point,(x0,y0).

Next value:

To get the next value $Y$ 2, we would use the value we just found for

Y

as follows:

Where
Y2 is the next estimated solution value;
Y1 is the current value;
$h$ is the interval between steps;
$X$ 1=X0+h; and
$f$ (x1,y1) is the value of the derivative at the current (x1,y1) point.

We continue this process for as many steps as required.

C Program for Euler's Method

#include<stdio.h>
#include<math.h>
#include<conio.h>
#define f(x,y) 2*(x)+(y)
//#define f(x,y) 3*(x)*(x)+1

int main(){
 float x,xp,x0,y0,y,h,m1,m2,m3,m4;
 printf("\n Euler's method'\n");
 printf("\nEnter initial values of x and y\n");
 scanf("%f%f",&x0,&y0);
 printf("\nEnter x at which function to be evaluated\n");
 scanf("%f",&xp);
 printf("\nEnter the step size \n");
 scanf("%f",&h);
 y=y0;
 x=x0;
 for(x=x0;x<xp;x=x+h){
  m1=f(x,y);
  printf("\nm1=%f",m1);
  y=y+h*(m1);
  printf("\ny=%f",y);
 }
 
 printf("Function value at x=%f is %f\n",xp,y);
 getch();
 return 0; 
 
}

Numerical Methods

Runge Kutta Method Formula and C Program

Arogya Thapa Magar at March 20, 2019 0

Runge Kutta Method

A Runge Kutta Method of numerically integrating ordinary differential equations by using a trial step at the midpoint of an interval to cancel out lower-order error terms. The second-order formula is

Each of these slope estimates can be described verbally.

k₁ is the slope at the beginning of the time step.
If we use the slope k₁ to step halfway through the time step, then k₂ is an estimate of the slope at the midpoint. This is the same as the slope, k₂, from the second-order midpoint method. This slope proved to be more accurate than k₁ for making new approximations for y(t).
If we use the slope k₂ to step halfway through the time step, then k₃ is another estimate of the slope at the midpoint.
Finally, we use the slope, k₃, to step all the way across the time step (to t₀+h), and k₄ is an estimate of the slope at the endpoint.

C Code For 4th Order Runge Kutta Method

#include<stdio.h>
#include<math.h>
#include<conio.h>
#define f(x,y) 1-2*(x)*(x)*(y)

int main(){
 float x,xp,x0,y0,y,h,m1,m2,m3,m4;
 printf("\n runge-Kutta Forth order\n");
 printf("\nEnter initial values of x and y\n");
 scanf("%f%f",&x0,&y0);
 printf("\nEnter x at which function to be evaluated\n");
 scanf("%f",&xp);
 printf("\nEnter the step size \n");
 scanf("%f",&h);
 y=y0;
 x=x0;
 for(x=x0;x<xp;x=x+h){
  m1=f(x,y);
  printf("\nm1=%f",m1);
  m2=f(x+h/2,y+(h*m1)/2);
  printf("\nm2=%f",m2);
  m3=f(x+h/2,y+(h*m2)/2);
  printf("\nm3=%f",m3);
  m4=f(x+h,y+h*m3);
  printf("\nm4=%f",m4);
  y=y+h/6*(m1+2*m2+2*m3+m4);
  printf("\ny=%f",y);
 }
 
 printf("Function value at x=%f is %f\n",xp,y);
 getch();
 return 0; 
 
}

C Program for Shooting Method | C Programming

Arogya Thapa Magar at March 19, 2019 0

Introduction

The shooting method is a method for solving a Boundary Value Problem by reducing it to the solution of an Initial Value Problem. Roughly speaking, we 'shoot' out trajectories in different directions until we find a trajectory that has the desired boundary value.

C Code For Shooting Method

#include<stdio.h>
#include<conio.h>
#include<math.h>
#define f1(x,y,z) (z)
#define f2(x,y,z) 6*(x)
int main(){
 float xa,xb,ya,yb,x0,y0,z0,x,y,z,xp,h,sol,ny,nz,error,E,g[3],v[3],gs;
 int i;
 printf("\nShooting method");
 printf("\nEnter the boundary conditions\n");
 scanf("%f%f%f%f",&xa,&ya,&xb,&yb);
 printf("Enter x at which value is required\n");
 scanf("%f",&xp);
 printf("Enter the step size\n");
 scanf("%f",&h);
 printf("Enter the accuracy limit\n");
 scanf("%f",&E);
 
 for(i=1; i<=2;i++){
  x=xa;
  y=ya;
  g[i]=z=i*(yb-ya)/(xb-xa);
  printf("g=%f\n",g[i]);
  while(x<xb){
   ny=y+(f1(x,y,z))*h;
   nz=z+(f2(x,y,z))*h;
   x=x+h;
   y=ny;
   z=nz;
   if(x==xp)
    sol=y;
  }
  v[i]=y;
  error=fabs(y-yb)/y;
  if(error<E){
   printf("y(%f)=%f",xp,sol);
   break;
  }
 }
 
 while(1){
  x=xa; y=ya;
  gs=z=g[2]-(v[2]-yb)/(v[2]-v[1])*(g[2]-g[1]);
  while(x<xb){
   ny=y+(f1(x,y,z))*h;
   nz=z+(f2(x,y,z))*h;
   x=x+h;
   y=ny;
   z=nz;
   if(x==xp)
    sol=y;
   
  }
  error=fabs(y-yb)/y;
  v[1]=v[2]; v[2]=y;
  g[1]=g[2]; g[2]=gs;
  
  if(error<E){
   printf("y(%f)=%f",xp,sol);
   break;
  }
 }
 
 getch();
 return 0;
}

Numerical Methods

C Program for Laplace Algorithm | C Programming

Arogya Thapa Magar at March 05, 2019 0

Laplace Algorithm

Laplace’s approximation is one of the most fundamental asymptotic techniques for the estimation of integrals containing a large parameter or variable. For integrals of the form

where f(t) and ψ(t) are real continuous functions defined on the interval [a,b](which may be infinite), the dominant contribution as x→+∞ arises from a neighborhood of the point where ψ(t) attains its maximum value. When ψ(t)possesses a single maximum at the point t0∈(a,b), so that ψ′(t0)=0, ψ″(t0)<0 span="">and f(t0)≠0, then I(x) has the asymptotic behavior

Abstract

A new algorithm for the numerical inversion of Laplace transforms is presented. The algorithm requires that transform values F(z) are available at arbitrary points in the complex z plane. It involves the approximation of the exponential function est in the inversion integral by a two-parameter function. The approximated inversion integral is then expanded into two different infinite series:

partial fraction expansion by the Mittag-Leffler theorem;
power series expansion in est.

In the numerical procedures, the summing of the infinite series are accelerated by a Wynn's type of in -algorithm. Trial runs on several test functions show that the results are highly accurate. The programming effort for the algorithm is minimal in comparison with other algorithms of comparable accuracy. The computation involves only a finite small number of function evaluations of the transform function. The method is believed to work for a wider range of functions since the algorithm does not make any assumption on representing f(t) by a linear combination of a certain class of basis functions.

C Code For Laplace Algorithm

#include<stdio.h>
#include<conio.h>
#include<math.h>
int main(){
	int n,i,j,k;
	float sum,error,E[10],a[10][10],b[10],new_x[10],old_x[10],tl,tr,tu,tb;
	printf("\nLaplace \n");
	printf("Enter the dimension of plate\n");
	scanf("%d",&n);
	printf("\nEnter the temperatures at left, right, bottom and top of the plate \n");
	scanf("%f%f%f%f",&tl,&tr,&tb,&tu);
	//construction of coefficient matrix
	for(i=0;i<=n;i++)
		a[i][i]=-4;
	for(i=0;i<=n;i++)
		a[i][n-i]=0;
	for(i=0;i<=n;i++)
	for(j=0;j<=n;j++){
		if(i!=j && j!=(n-i))
			a[i][j]=1;
	}
	
	//Construction of RHS vector
	for(i=0;i<=n;i++)
		b[i]=0;
		
	k=0;
	for(i=1;i<n;i++){
		for(j=1;j<n;j++){
			if((i-1)==0)
				b[k]=b[k]-tl;
			if((i+1)==n)
				b[k]=b[k]-tr;
			if((j-1)==0)
				b[k]=b[k]-tb;
			if((j+1)==n)
				b[k]=b[k]-tu;
			k++;
		}
	}
	
	printf("\nEnter accuracy limit\n");
	scanf("%f",&error);
	
	//Solving the system of linear equations using Gauss-Seidel method
	for(i=0;i<=n;i++){
		new_x[i]=0;
	}
	
	while(1){
		for(i=0;i<=n;i++){
			sum=b[i];
			for(j=0;j<=n;j++){
				if(i!=j)
					sum=sum-a[i][j]*new_x[j];
			}
			old_x[i]=new_x[i];
			new_x[i]=sum/a[i][i];
			E[i]=fabs(new_x[i]-old_x[i])/fabs(new_x[i]);
		}
		for(i=0;i<=n;i++){
			if(E[i]>error)
				break;
				
		}
		if(i==(n+1))
			break;
		else
			continue;
			
	}
	
	printf("\nsolution\n");
	for(i=0;i<=n;i++)
	printf("x%d=%6.2f\n",i+1,new_x[i]);
	getch();
	return 0;
}

C Program To Find The Value Using Poison Algorithm | C Programming

Arogya Thapa Magar at March 04, 2019 0

C Programming | C Program to find the value using poison algorithm

#include<stdio.h>
#include<conio.h>
#include<math.h>
#define g(x,y) 2*(x)*(x)*(y)*(y)
int main(){
 int n,i,j,k;
 float sum,error,E[10],a[10][10],b[10],new_x[10],old_x[10],tl,tr,tu,tb,h,buffer;
 printf("\nPoision's' \n");
 printf("Enter the dimension of plate\n");
 scanf("%d",&n);
 printf("\nEnter the dimension of grid\n");
 scanf("%f",&h);
 printf("\nEnter the temperatures at left, right, bottom and top of the plate \n");
 scanf("%f%f%f%f",&tl,&tr,&tb,&tu);
 //construction of coefficient matrix
 for(i=0;i<=n;i++)
  a[i][i]=-4;
 for(i=0;i<=n;i++)
  a[i][n-i]=0;
 for(i=0;i<=n;i++)
 for(j=0;j<=n;j++){
  if(i!=j && j!=(n-i))
   a[i][j]=1;
 }
 
 //Construction of RHS vector
 k=0;
 for(i=1;i<n;i++)
 for(j=1;j<n;j++)
  b[k++]=h*h*g(i,j);
 
 k=0;
 for(i=1;i<n;i++){
  for(j=1;j<n;j++){
  
   if((i-1)==0)
    b[k]=b[k]-tl;
    
   if((i+1)==n)
    b[k]=b[k]-tr;
    
   if((j-1)==0)
    b[k]=b[k]-tb;
   
   if((j+1)==n)
    b[k]=b[k]-tu;
    
   k++; 
  }
  
 }
 
 printf("\nEnter accuracy limit\n");
 scanf("%f",&error);
 
 //Solving the system of linear equations using Gauss-Seidel method
 for(i=0;i<=n;i++){
  new_x[i]=0;
 }
 
 while(1){
  for(i=0;i<=n;i++){
   sum=b[i];
   for(j=0;j<=n;j++){
    if(i!=j)
     sum=sum-a[i][j]*new_x[j];
   }
   old_x[i]=new_x[i];
   new_x[i]=sum/a[i][i];
   E[i]=fabs(new_x[i]-old_x[i])/fabs(new_x[i]);
  }
  for(i=0;i<=n;i++){
   if(E[i]>error)
    break;
    
  }
  if(i==(n+1))
   break;
  else
   continue;
   
 }
 
 printf("\nsolution\n");
 for(i=0;i<=n;i++)
 printf("x%d=%6.2f\n",i+1,new_x[i]);
 getch();
 return 0;
}