Solving simultaneous equations

Determinant

1. 2 ×2 matrices
   

   ⎢ ⎢ ⎢ ⎢ a11 a12 a21 a22. ⎢ ⎢ ⎢ ⎢ ≡ a11 a22− a12 a21

   (4)
2. 3 ×3 matrices
   

   ⎢ ⎢ ⎢ ⎢ ⎢ a11 a12 a13 a21 a22 a23 a31 a32 a33 ⎢ ⎢ ⎢ ⎢ ⎢ ≡ a11 a22 a33 + a12 a23 a31+ a21 a32 a13− a13 a22a31−a12 a21 a33 − a23 a32 a11

   (5)
3. The determinant of n ×n (n > 3) matrices can be expressed similar to the above. It consists of n! terms each of which is a product of n variables.

Some properties of determinants

1. | A B | = | A | | B |
2. | I | = 1 where I is the identity matrix.
3. | A |=0 ⇒ A is singular.

Determinant/Cramer's rule

Gauss-Jordan elimination method

Inverse matrix by the Gauss-Jordan elimination method (optional)

Implementation

#include <stdio.h>
#define N 3
int main()
{
 double a[N][N+1]={{2, 3, 4, 6},{3, 5, 2, 5},{4, 3, 30, 32}};
  double pivot,d;
  int i,j,k;

for(k=0; k<N; k++)
{
 pivot=a[k][k];

 for(j=k; j<N+1; j++) a[k][j]=a[k][j]/pivot;
 for(i=0; i<N;  i++)
  {
    if(i != k)
     {
       d=a[i][k]; 
        for(j=k; j<N+1; j++) a[i][j]=a[i][j]-d*a[k][j];
     }
   }
}

for(i=0; i<N; i++) printf("x[%d]=%lf\n", i+1, a[i][N]);
return 0;
}

LU decomposition (optional)

1. Solve L y = c first.
2. Then, solve U x = y.

Notes

1. If A is a symmetric matrix (a_ij=a_ji), LU-decomposition is called the Cholesky decomposition. The number of operations required for the Cholesky decomposition is n³/6.
2. If A is a sparse matrix, so are L and U.
3. If A is a triangular diagonal matrix (typical of the finite element method), so are L and U.

Gauss-Seidel method (Important)

#include <stdio.h>
int main()
{
double x, y, z;
int i,n;

x=y=z=0.0;

printf("Enter # of iteration = ");
scanf("%d", &n);

for (i=0;i<n;i++)
{
x = (10-y-2*z)/7;
y = (8-x-3*z)/8.0;
z = (6-2*x-3*y)/9.0;
}

printf("x = %lf, y= %lf, z=%lf\n", x,y,z);
return 0;
}

Numerical differentiation

Forward difference

Backward difference

Central difference

• Forward difference:
  

  f ′(x) ≈ f ( x + h )−f (x) h + O(h).
• Backward difference:
  

  f ′(x) ≈ f (x)−f ( x − h ) h + O(h).
• Central difference:
  

  f ′(x) ≈ f ( x + h )−f ( x − h ) 2h + O(h2).

Second order derivative

Example

#include <stdio.h>
#define N 11

int main()
{
double y[N]={0, 0.0998, 0.1986, 0.2955, 0.3894, 0.4794, 0.5646, 0.6442, 0.7173, 0.7833, 0.8414};
double central[N], h=0.1;
int i;

/* Excluding x=0 and x=1.0 */
for (i=1;i<N-1;i++) central[i]= (y[i+1]-y[i-1])/(2*h);

printf ("   x    Central \n---------------------------\n");

for (i=1;i<N-1;i++)
 printf ("%lf %lf\n", i*h, central[i]);

return 0;
}

What if you want to obtain f′(x) at x = a with the same accuracy as the central difference ?

#include <stdio.h>
#define N 11

int main()
{
double y[N]={0, 0.0998, 0.1986, 0.2955, 0.3894, 0.4794, 0.5646, 0.6442, 0.7173, 0.7833, 0.8414};
double central[N], h=0.1;
int i;

for (i=1;i<N-1;i++) central[i]= (y[i+1]-y[i-1])/(2*h);

central[0]=(4*y[1]-y[2]-3*y[0])/(2*h);

printf ("   x    Central \n---------------------------\n");

for (i=0;i<N-1;i++)
 printf ("%lf %lf\n", i*h, central[i]);

return 0;
}

Footnotes: