Numerical Method

bisection

bisection[f_, l_ : 0, u_ : 1, iters_ : 10] := (
  a = l;
  b = u;
  If[f[a]*f[b] > 0,
   Print["We cannot find the roots in the given interval"],
   
   c = N[(a + b) / 2];
   OutputDetails = {{0, a, b, c, f[c]}};
   
   Do[
    If[f[a]*f[c] > 0, a = c, b = c];
    c = N[(a + b)/2];
    OutputDetails = Append[OutputDetails, {i, a, b, c, f[c]}];
    , {i, iters}];
   
   Print[NumberForm[TableForm[OutputDetails,
      TableHeadings -> {None, {"i", "ai", "bi", "ci", "f[ci]"}}], 
     6]];
   
   Print["Root is ", NumberForm[c, 6]];
   
   Plot[f[x], {x, l, u}, PlotStyle -> {Red}, 
    AxesLabel -> {"x", "f(x)"}]
   ]
  )

euler

euler[eqn_, point_, target_, iters_ : 10] := (
  {xn, yn} = point;
  dydx = Values[Solve[eqn, y'[x]][[1]][[1]]];
  h = N[(target - xn)/iters];
  OutputDetails = {{0, xn, yn}};
  Do[
   yn = yn + h*(dydx /. {x -> xn, y[x] -> yn});
   xn = xn + h;
   OutputDetails = Append[OutputDetails, {i, xn, yn}];
   , {i, iters}
   ];
  Print[NumberForm[TableForm[OutputDetails,
     TableHeadings -> {None, {"i", "xi", "yi"}}], 6]];
  Print["Value of y[", target, "] = ", yn];
  )

gauss-elimination

gaussianElimination[eqns_, vars_] := (
  {b, A} = Normal@CoefficientArrays[eqns, vars];
  Print["A = ", MatrixForm[A], ", b = ", MatrixForm[b]];
  Ab = ArrayFlatten[{{A, Transpose[{b}]}}];
  Print["Augmented Matrix : ", MatrixForm[Ab]];
  Do[
   For[j = i + 1, j <= Length[Ab], j++,
    Print["R", j, " -> R", j, " - ", A[[j]][[i]], " x R", i];
    Ab[[j]] = Ab[[j]] - (Ab[[j]][[i]]/Ab[[i]][[i]]) Ab[[i]];
    ];
   Print["\t\t\t\t=> ", MatrixForm[Ab]];
   , {i, Length[Ab] - 1}
   ];
  Do[
   Ab[[i]] = Ab[[i]]/Ab[[i]][[i]];
   Do[
    If[i != j,
     Ab[[j]] = Ab[[j]] - Ab[[j]][[i]] Ab[[i]],
     Null]
    , {j, Length[Ab]}
    ];
   , {i, Length[Ab]}
   ];
  Print[vars == Transpose[Ab][[-1]]*-1];
  )

gauss-jordan

gaussianJordan[eqns_, vars_] := (
  {b, A} = Normal@CoefficientArrays[eqns, vars];
  Print["A = ", MatrixForm[A], ", b = ", MatrixForm[b]];
  Ab = ArrayFlatten[{{A, Transpose[{b}]}}];
  Print["Augmented Matrix : ", MatrixForm[Ab]];
  Do[
   Print["R", i, " -> R", i, " / ", Ab[[i]][[i]]];
   Ab[[i]] = Ab[[i]]/Ab[[i]][[i]];
   Do[
    If[i != j,
     Print["R", j, " -> R", j, " - ", A[[j]][[i]], " x R", i];
     Ab[[j]] = Ab[[j]] - Ab[[j]][[i]] Ab[[i]],
     Null]
    , {j, Length[Ab]}
    ];
   Print["\t\t\t\t=> ", MatrixForm[Ab]];
   , {i, Length[Ab]}
   ];
  Print[vars == Transpose[Ab][[-1]]*-1];
  )

jacobi

jacobi[eqns_, vars_, iters_ : 10, errbound_ : 0, 
  iv_ : ConstantArray[0, 100]] := (
  xs = {};
  xvs = {};
  Do[
   xs = Append[xs, Solve[eqns[[i]], vars[[i]]][[1]][[1]]];
   xvs = Append[xvs, vars[[i]] -> iv[[i]]];
   , {i, Length[vars]}
   ];
  OutputDetails = {Join[{0}, Values[xvs]]};
  maxnorm = Infinity;
  For[i = 1, maxnorm >= errbound && i <= iters, i++,
   txvs = xvs;
   Do[
    xvs[[j]] = 
      xs[[j]] /. Select[txvs, Position[txvs, #][[1]][[1]] != j &] ;
    , {j, Length[vars]}
    ];
   maxnorm = Max[Abs[Values[xvs] - Values[txvs]]];
   OutputDetails = Append[OutputDetails, Join[{i}, N[Values[xvs]]]];
   ];
  Print[NumberForm[TableForm[OutputDetails,
     TableHeadings -> {None, Join[{"i"}, vars]}], 6]];
  Print[N[xvs]];
  )

lagrange

lagrangeInterpolation[xs_, ys_] := (
  Assert[Length[xs] == Length[ys]];
  n = Length[xs];
  V  = {Table[x^i, {i, n, 0, -1}]};
  Do[
   V = Append[V, Join[xs[[i]]^Range[n, 1, -1], {1}]];
   , {i, n}
   ];
  V[[1]][[1]] = y[x];
  Do[
   V[[i]][[1]] = ys[[i - 1]]
   , {i, 2, n + 1}
   ];
  
  Print["Vandermonde's Determinant: ", MatrixForm[V]];
  yy = Solve[Det[V] == 0, y[x]];
  Print["Interpolated polynomial : ", yy[[1]][[1]]];
  Show[Plot[Values[yy[[1]][[1]]], {x, Min[xs] - 3, Max[xs] + 3}],
   ListPlot[Transpose@{xs, ys}, 
    PlotStyle -> {Red, PointSize[Large]}]]
  )

newton-interpolation

newtonDividedDifference[xs_, ys_] := (
  dely = {ys};
  delx = {xs};
  Do[
   dy = dely[[-1]];
   dx = delx[[-1]];
   dy = dy[[2 ;;]] - dy[[;; -2]];
   dx = dx[[2 ;;]] - dx[[;; -2]];
   dely = Append[dely, dy];
   delx = Append[delx, dx];
   , {i, Length[xs] - 1}
   ];
  yy = 0;
  Do[
   yy = yy + Product[(x - xs[[j]]), {j, 1, i - 1}]* dely[[i]][[1]];
   , {i, Length[xs]}
   ];
  Print[y[x] == yy];
  Show[Plot[yy, {x, Min[xs] - 3, Max[xs] + 3}],
   ListPlot[Transpose@{xs, ys}, 
    PlotStyle -> {Red, PointSize[Large]}]]
  )

newtonRaphson

newtonRaphson[f_, iters_ : 10, init_ : 1, mnfx_ : 0] := (
  x = init;
  OutputDetails = {{0, x, f[x]}};
  
  For[i = 1, Abs[f[x]] >= mnfx && i <= iters, i++,
   x = x - N[f[x]/f'[x]];
   OutputDetails = Append[OutputDetails, {i, x, f[x]}];
   ];
  
  Print[NumberForm[TableForm[OutputDetails,
     TableHeadings -> {None, {"i", "x", "f[x]"}}], 6]];
  
  Print["Root is ", NumberForm[x, 6]];
  
  Plot[f[xx], {xx, x - 10, x + 10} , PlotStyle -> {Red}, 
   AxesLabel -> {"x", "f(x)"}]
  )

regulaFalsi

regulaFalsi[f_, l_, u_, iters_ : 10, errbound_ : 0 ] := (
  {a, b} = {l, u};
  
  If [f[a]*f[b] > 0,
   Print["We cannot find the roots in the given interval"],
   
   c = N[(a*f[b] - b*f[a])/(f[b] - f[a])];
   OutputDetails = {{0, a, b, c, f[c]}};
   
   For[i = 1, Abs[f[c]] >= errbound && i <= iters, i++,
    If[f[a]*f[c] < 0, b = c, a = c];
    c = N[(a*f[b] - b*f[a])/(f[b] - f[a])];
    OutputDetails = Append[OutputDetails, {i, a, b, c, f[c]}];
    ];
   
   Print[NumberForm[TableForm[OutputDetails,
      TableHeadings -> {None, {"i", "ai", "bi", "ci", "f[ci]"}}], 
     6]];
   
   Print["Root is ", NumberForm[c, 6]];
   
   Plot[f[x], {x, l, u} , PlotStyle -> {Red}, 
    AxesLabel -> {"x", "f(x)"}]
   ]
  )

romberg

romberg[f_, lowl_, upl_, m_] := (
  OutputDetails = {};
  Do[ 
   h = N[(upl - lowl)/Power[2, n - 1]];
   ys = f[lowl + Range[1, Power[2, n - 1] - 1]*h];
   t[n, 1] = (h/2)*(f[lowl] + f[upl] + 2*Total[ys]);
   OutputDetails = Append[OutputDetails, {n, h, N[t[n, 1], 12]}];
   , {n, m}];
  Print[NumberForm[TableForm[OutputDetails,
     TableHeadings -> {None, {"n", "h", "Trapezoidal approx"}}], 6]];
  ti = N[Integrate[f[x], {x, 0, 1}], 12];
  Print["Actual value of the integral is: ", ti];
  OutputDetails = {};
  Do[ 
    Do[ 
     t[i, j] = 
      N[(Power[4, j - 1]*t[i, j - 1] - 
          t[i - 1, j - 1])/(Power[4, j - 1] - 1), 12];
     OutputDetails = Append[OutputDetails, {i, j, t[i, j]}]
     , {j, 2, i}]
    , {i, 2, m}]
   Print[NumberForm[TableForm[OutputDetails,
      TableHeadings -> {None, {"i", "j", "Romberg approx, t[i,j]"}}], 
     6]];
  Print["The Value of the Integral by Romberg Integration is: ", 
   t[m, m]];
  Print["Error is: ", Abs[ti - t[m, m]]]
  )

secant

secant[f_, a_ : 0, b_ : 1, iters_ : 10, errbound_ : 0] := (
  {xp, xn} = {a, b};
  OutputDetails = {{0, xp, xn, f[xn]}};
  For[i = 0, Abs[f[xn]] >= errbound && i <= iters, i++,
   {xp, xn} = {xn, N[(xp*f[xn] - xn*f[xp])/(f[xn] - f[xp])]};
   OutputDetails = Append[OutputDetails, {i, xp, xn, f[xn]}];
   ];
  
  Print[NumberForm[TableForm[OutputDetails,
     TableHeadings -> {None, {"i", "x[n-1]", "x[n]", "f[x[n]]"}}], 6]];
  
  Print["Root is ", NumberForm[xn, 6]];
  
  Plot[f[x], {x, a, b}, PlotStyle -> {Red}, AxesLabel -> {"x", "f(x)"}]
  )

seidel

gaussSeidel[eqns_, vars_, iters_ : 10, errbound_ : 0, 
  iv_ : ConstantArray[0, 100]] := (
  xs = {};
  xvs = {};
  Do[
   xs = Append[xs, Solve[eqns[[i]], vars[[i]]][[1]][[1]]];
   xvs = Append[xvs, vars[[i]] -> iv[[i]]];
   , {i, Length[vars]}
   ];
  OutputDetails = {Join[{0}, Values[xvs]]};
  maxnorm = Infinity;
  For[i = 1, maxnorm >= errbound && i <= iters, i++,
   maxnorm = 0;
   Do[
    nv = xs[[j]] /. Select[xvs, Position[xvs, #][[1]][[1]] != j &] ;
    maxnorm = Max[{maxnorm, Abs[Values[xvs[[j]]] - Values[nv]]}];
    xvs[[j]] = nv;
    , {j, Length[vars]}
    ];
   OutputDetails = Append[OutputDetails, Join[{i}, N[Values[xvs]]]];
   ];
  Print[NumberForm[TableForm[OutputDetails,
     TableHeadings -> {None, Join[{"i"}, vars]}], 6]];
  Print[N[xvs]];
  )

simpson

simpson[f_, lowl_, upl_, n_] := (
  h = (upl - lowl)/n;
  ys = f[Range[lowl, upl, h]];
  i = (h/3) (ys[[1]] + ys[[-1]] + 
      2*Total[ys[[3 ;; Length[ys] - 1 ;; 2]]] + 
      4*Total[ys[[2 ;; Length[ys] - 1 ;; 2]]]);
  
  Print["Simpson integral = ", TraditionalForm[i], " = ", N[i]];
  
  ti = Integrate[f[x], {x, lowl, upl}];
  Print["True Integral, Itrue = ", TraditionalForm[ti], " = ", 
   N[ti]];
  Print["Absolute error for simpson rule = |Simpson Approx - True \
Value| = ", Abs[N[i - ti]]]
  )

trapezoidal

trapezoidal[f_, lowl_, upl_, n_] := (
  h = (upl - lowl)/n;
  ys = f[Range[lowl, upl, h]];
  i = h*((ys[[1]] + ys[[-1]])/2 + Total[ys[[2 ;; Length[ys] - 1]]]);
  Print["Approx Integral, Iapprox = ", TraditionalForm[i], " = ", 
   N[i]];
  
  ti = Integrate[f[x], {x, lowl, upl}];
  Print["True Integral, Itrue = ", TraditionalForm[ti], " = ", 
   N[ti]];
  Print["Absolute error = |Trapezoidal Approx - True Value| = ", 
   Abs[N[i - ti]]]
  )