Bisection method
Bisection[a0_,b0_,e0_,n_]:=
Module[{},
a=N[a0];
b=N[b0];
e=N[e0];
i=0;
output={{a,b,i,}};
While[i<n,
c=(a+b)/2;
output=Append[output,{i+1,a,b,c}];
If[Sign[f[b]]==Sign[f[c]],b=c,a=c];
If[(b-a)/2<e,Print["Condition exists at ",i+1," . "];Break[]];
i=i+1;
];
Print[NumberForm[TableForm[output,TableHeadings->{None,{"i","a{i}","b{i}","c{i}"}}],16]];
Print[" Root p = ", NumberForm[c,16]];
]
Newton method
NewtonMethod[p0_,e0_,n_]:=
Module[{},
p=N[p0];
e=N[e0];
i=0;
output={{i,p}};
While[i<n,
p1=p-f[p]/f'[p];
output=Append[output,{i+1,p1}];
If[Abs[p1-p]<e,Print["Condition exists at ",i+1," . "];Break[]];
p=p1;
i=i+1;
];
Print[NumberForm[TableForm[output,{"i","p"}],18]];
Print[" p = ",NumberForm[p1,18]];
]
Secant method
SecantMethod[a0_,b0_,e0_,n_]:=
Module[{},
p=N[a0];
p1=N[b0];
e=N[e0];
i=0;
output={{i,p}};
While[i<n,
p2=p1-f[p1]*(p1-p)/(f[p1]-f[p]);
output=Append[output,{i+1,p2}];
If[Abs[p2-p1]<e,Print["Condition exists at ",i+1," . "];Break[]];
p=p1;
p=p2;
i=i+1;
];
Print[NumberForm[TableForm[output,{"i","p"}],18]];
Print[" p = ",NumberForm[p2,18]];
]
LUDecomposition(DooLittle and Crouts)
LUDoolittle[Data_]:=Module[{},
ClearAll;
U={{0,0,0},{0,0,0},{0,0,0}};
L={{1,0,0},{0,1,0},{0,0,1}};
n=Length[Data];
For[i=1,i<=n,i++,U[[1,i]]=Data[[1,i]]];
For[j=2,j<=n,j++,L[[j,1]]=Data[[j,1]]/U[[1,1]]];
For[k=2,k<=n-1,k++,
For[i=k,i<=n,i++,
U[[k,i]]=Data[[k,i]]-Sum[L[[k,j]]*U[[j,i]],{j,1,k-1}]
]
]
For[k=2,k<=n,k++,
For[j=k+1,j<=n,j++,
L[[j,k]]=(Data[[j,k]]-Sum[L[[j,i]]*U[[i,k]],{i,1,k-1}])/U[[k,k]]
]
]
For[k=3,k<=n,k++,
For[i=k,i<=n,i++,
U[[k,i]]=Data[[k,i]]-Sum[L[[k,j]]*U[[j,i]],{j,1,k-1}]
]
]
Print["L = ",MatrixForm[L]];
Print["U = ",MatrixForm[U]];
Print["Result Verification : ",MatrixForm[L],MatrixForm[U]," = ",MatrixForm[L.U]];
]
LUCrout[Data_]:=Module[{},
ClearAll;
L={{0,0,0},{0,0,0},{0,0,0}};
U={{1,0,0},{0,1,0},{0,0,1}};
n=Length[Data];
For[i=1,i<=n,i++,L[[i,1]]=Data[[i,1]]];
For[j=2,j<=n,j++,U[[1,j]]=Data[[1,j]]/L[[1,1]]];
For[k=2,k<=n-1,k++,
For[i=k,i<=n,i++,
L[[i,k]]=Data[[i,k]]-Sum[L[[i,j]]*U[[j,k]],{j,1,k-1}]
]
]
For[k=2,k<=n,k++,
For[j=k+1,j<=n,j++,
U[[k,j]]=(Data[[k,j]]-Sum[L[[k,i]]*U[[i,j]],{i,1,k-1}])/L[[k,k]]
]
]
For[k=3,k<=n,k++,
For[i=k,i<=n,i++,
L[[i,k]]=Data[[i,k]]-Sum[L[[i,j]]*U[[j,k]],{j,1,k-1}]
]
]
Print["L = ",MatrixForm[L]];
Print["U = ",MatrixForm[U]];
Print["Result Verification : ",MatrixForm[L],MatrixForm[U]," = ",MatrixForm[L.U]];
]
Gauss-Jacobi method
GaussJacobi[a0_,b0_,x0_,e0_,m0_]:=
Module[{},
A=a0;
B=b0;
X=x0;
n=Length[X];
e=N[e0];
m=N[m0];
k=0;
X1=X;
Print["Given System is : ",MatrixForm[A],"X = ",MatrixForm[B]];
output={{k,NumberForm[X1,10]}};
While[k<m,
For[i=1,i<=n,i++,
X1[[i]]=N[(1/A[[i,i]])*(B[[i]]-Sum[A[[i,j]]*X[[j]],{j,1,i-1}]-Sum[A[[i,j]]*X[[j]],{j,i+1,n}])];
];
output=Append[output,{k+1,NumberForm[X1,10]}];
If[Norm[X1-X]<e,Print["Condition Exists at ",k+1," . "];Break[]];
X=X1;
k++
];
Print[NumberForm[TableForm[output, TableHeadings->{None,{"k","X[]"}}],16]];
Print["X = ",NumberForm[X1,10]];
]
Gauss-Seidel method
GaussSeidel[a0_,b0_,x0_,e0_,m0_]:=
Module[{},
A=a0;
B=b0;
X=x0;
n=Length[X];
e=N[e0];
m=N[m0];
k=0;
X1=X;
Print["Given System is : ",MatrixForm[A],"X = ",MatrixForm[B]];
output={{k,NumberForm[X1,10]}};
While[k<m,
For[i=1,i<=n,i++,
X1[[i]]=N[(1/A[[i,i]])*(B[[i]]-Sum[A[[i,j]]*X1[[j]],{j,1,i-1}]-Sum[A[[i,j]]*X[[j]],{j,i+1,n}])];
];
output=Append[output,{k+1,NumberForm[X1,10]}];
If[Norm[X1-X]<e,Print["Condition Exists at ",k+1," . "];Break[]];
X=X1;
k++
];
Print[NumberForm[TableForm[output,TableHeadings->{None,{"k","X[]"}}],16]];
Print["X = ",NumberForm[X1,10]];
]
Lagrange and NDD
Lagrange[Data_]:=Module[{},
X[k_]:=Transpose[Data][[1,k+1]];
Y[k_]:=Transpose[Data][[2,k+1]];
n=Length[Data]-1;
For[k=0,k<=n,k++,
L[n,k,x_]:=Expand[Product[(x-X[j])/(X[k]-X[j]),{j,0,k-1}]*Product[(x-X[j])/(X[k]-X[j]),{j,k+1,n}]]
];
Return[Expand[Sum[Y[k]*L[n,k,x],{k,0,n}]]];
]
NDD[Data_,i_,j_]:=Module[{},
X[k_]:=Transpose[Data][[1,k+1]];
Y[k_]:=Transpose[Data][[2,k+1]];
If[i==j,Return[Y[i]],
Return[(NDD[Data,i+1,j]-NDD[Data,i,j-1])/(X[j]-X[i])]];
];
NewtonInterpolation[Data_]:=Module[{},
X[k_]:=Transpose[Data][[1,k+1]];
Y[k_]:=Transpose[Data][[2,k+1]];
n=Length[Data]-1;
Return[Expand[Sum[NDD[Data,0,k]*Product[(x-X[i]),{i,0,k-1}],{k,0,n}]]];
];
Trapezoidal Rule Code
a=1;
b=2;
h=(b-a)/2;
f[x_]:=x^5+2x^4+x+1;
SI=(h)*(f[a]+f[b]);
Print["Integration by Trapezoidal rule : ",N[SI]];
DI=Integrate[f[x],{x,1,2}];
Print["Integration by direct rule : ",N[DI]];
Print["Error : ",Abs[N[SI]-N[DI]]]
Simpson’s 1/3 rule
a=1;
b=2;
h=(b-a)/2;
f[x_]:=Exp[x^2]+Exp[-x];
SI=(h/3)*(f[a]+4f[a+h]+f[a+2h]);
Print["Integration by Simpson's 1/3rd rule : ",N[SI]];
DI=Integrate[f[x],{x,1,2}];
Print["Integration by direct rule : ",N[DI]];
Print["Error : ",Abs[N[SI]-N[DI]]]
Simpson’s 3/8 rule
a=1;
b=2;
h=(b-a)/3;
f[x_]:=x^5+2x^4+x+1;
SI=(3h/8)*(f[a]+3f[a+h]+3f[a+2h]+f[a+3h]);
Print["Integration by Simpson's 3/8th rule : ",N[SI]];
DI=Integrate[f[x],{x,1,2}];
Print["Integration by direct rule : ",N[DI]];
Print["Error : ",Abs[N[SI]-N[DI]]]
Euler’s method
Eulermethod[a0_,b0_,h0_,f_,alpha_]:=Module[{a=N[a0],b=N[b0],h=N[h0],n,x},
n=(b-a)/h;
y[0]=alpha;
For[i=0,i<=n,i++,
x[i]=a+h*i;
y[i+1]=y[i]+h*f[x[i],y[i]];
Print["value at x[",i,"] = ",x[i]," is ",y[i]];
];
]
Runge-Kutta Code
RungeKutta[a0_,b0_,n0_,f_,alpha_]:=Module[{a=N[a0],b=N[b0],n=N[n0],h,ti,k1},
h=(b-a)/n;
ti=Table[a+(j-1)h,{j,1,n+1}];
wi=Table[0,{n+1}];
wi[[1]]=alpha;
outputdetails={{0,ti[[i]],alpha}};
For[i=1,i<=n,i++,
k1=wi[[i]]+(h/2)*f[ti[[i]],wi[[i]]];
wi[[i+1]]=wi[[i]]+h*f[ti[[i]]+h/2,k1];
outputdetails=Append[outputdetails,{i,N[ti[[i+1]]],N[wi[[i+1]]]}];];
Print[NumberForm[TableForm[outputdetails,TableHeadings->{None,{"i","ti","wi"}}],6]];]
