summaryrefslogtreecommitdiff
path: root/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao
diff options
context:
space:
mode:
authorPrashant S2020-04-14 10:25:32 +0530
committerGitHub2020-04-14 10:25:32 +0530
commit06b09e7d29d252fb2f5a056eeb8bd1264ff6a333 (patch)
tree2b1df110e24ff0174830d7f825f43ff1c134d1af /Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao
parentabb52650288b08a680335531742a7126ad0fb846 (diff)
parent476705d693c7122d34f9b049fa79b935405c9b49 (diff)
downloadall-scilab-tbc-books-ipynb-master.tar.gz
all-scilab-tbc-books-ipynb-master.tar.bz2
all-scilab-tbc-books-ipynb-master.zip
Merge pull request #1 from prashantsinalkar/masterHEADmaster
Initial commit
Diffstat (limited to 'Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao')
-rw-r--r--Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/1-Basics_in_Computing.ipynb119
-rw-r--r--Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/10-Elliptical_Partial_Differential_Equations.ipynb299
-rw-r--r--Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/11-Hyperbolic_Partial_Differential_Equations.ipynb137
-rw-r--r--Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/2-Solution_of_Algebraic_and_Transcendental_Equations.ipynb642
-rw-r--r--Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/3-Solution_of_Linear_System_of_Equations_and_Matrix_Inversion.ipynb549
-rw-r--r--Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/4-Eigenvalue_Problems.ipynb278
-rw-r--r--Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/5-Curve_Fitting.ipynb441
-rw-r--r--Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/6-Interpolation.ipynb1135
-rw-r--r--Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/7-Numerical_Differentiation_and_Integration.ipynb887
-rw-r--r--Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/8-Ordinary_Differential_Equations.ipynb497
-rw-r--r--Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/9-Parabolic_Partial_Differential_Equations.ipynb429
11 files changed, 5413 insertions, 0 deletions
diff --git a/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/1-Basics_in_Computing.ipynb b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/1-Basics_in_Computing.ipynb
new file mode 100644
index 0000000..cdff1d6
--- /dev/null
+++ b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/1-Basics_in_Computing.ipynb
@@ -0,0 +1,119 @@
+{
+"cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Chapter 1: Basics in Computing"
+ ]
+ },
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 1.1: Conversion_of_Decimal_to_Binary.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 1.1\n",
+"clc\n",
+"clear\n",
+"\n",
+"dec_N = 47;\n",
+"bin_N = dec2bin(dec_N)\n",
+"disp(bin_N)"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 1.2: Conversion_of_Binary_to_Decimal.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 1.2\n",
+"clc\n",
+"clear\n",
+"\n",
+"dec = 0.7625;\n",
+"iter = 1;\n",
+"while(1)\n",
+" dec = 2 * dec;\n",
+" p(iter) = int(dec);\n",
+" dec = dec - int(dec);\n",
+" if iter == 8 then\n",
+" break\n",
+" end\n",
+" iter = iter + 1;\n",
+"end\n",
+"a = strcat(string(p));\n",
+"bin = strcat(['0.',a])\n",
+"disp(bin)"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 1.3: Conversion_of_Decimal_to_Binary_and_Octal.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 1.3\n",
+"clc\n",
+"clear\n",
+"\n",
+"dec_N = 59;\n",
+"bin_N = dec2bin(dec_N)\n",
+"oct_N = dec2oct(dec_N)\n",
+"disp(bin_N,'Binary:')\n",
+"disp(oct_N,'Octal:')"
+ ]
+ }
+],
+"metadata": {
+ "kernelspec": {
+ "display_name": "Scilab",
+ "language": "scilab",
+ "name": "scilab"
+ },
+ "language_info": {
+ "file_extension": ".sce",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "scilab",
+ "version": "0.7.1"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/10-Elliptical_Partial_Differential_Equations.ipynb b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/10-Elliptical_Partial_Differential_Equations.ipynb
new file mode 100644
index 0000000..2c0a01f
--- /dev/null
+++ b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/10-Elliptical_Partial_Differential_Equations.ipynb
@@ -0,0 +1,299 @@
+{
+"cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Chapter 10: Elliptical Partial Differential Equations"
+ ]
+ },
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 10.1: Laplace_Equation_using_Five_Point_Formulae.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 10.1\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"h = 1/4;\n",
+"xf = 1;\n",
+"yf = 1;\n",
+"x = 0:h:xf;\n",
+"y = 0:h:yf;\n",
+"m = length(y);\n",
+"n = length(x);\n",
+"\n",
+"u = zeros(m,n);\n",
+"u(m,:) = 100*x;\n",
+"u(:,n) = 100*y';\n",
+"u0 = u;\n",
+"\n",
+"I = ceil(m/2);\n",
+"J = ceil(n/2);\n",
+"\n",
+"u(J,I) = (u0(J-2,I-2) + u0(J-2,I+2) + u0(J+2,I-2) + u0(J+2,I+2)) / 4;\n",
+"\n",
+"for j = [J-1 J+1]\n",
+" for i = [I-1 I+1]\n",
+" u(j,i) = (u(j-1,i-1) + u(j-1,i+1) + u(j+1,i-1) + u(j+1,i+1)) / 4;\n",
+" end\n",
+"end\n",
+"\n",
+"j1 = [J-1 J J J+1];\n",
+"i1 = [I I-1 I+1 I];\n",
+"for k = 1:4\n",
+" i = i1(k);\n",
+" j = j1(k);\n",
+" u(j,i) = (u(j,i-1) + u(j,i+1) + u(j-1,i) + u(j+1,i)) / 4;\n",
+"end\n",
+"\n",
+"disp(u,'u:')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 10.2: Temperature_in_Two_Dimensional_Geometry.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"// Example 10.2\n",
+"// This is an analytical problem and need not be coded."
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 10.3: Laplace_Equation_in_Two_Dimension_using_Five_Point_Formulae.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 10.3\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"m = 5;\n",
+"n = 5;\n",
+"u = zeros(m,n);\n",
+"u(m,:) = [50 100 100 100 50];\n",
+"u0 = u;\n",
+"I = ceil(m/2);\n",
+"J = ceil(n/2);\n",
+"\n",
+"u(J,I) = (u0(J-2,I-2) + u0(J-2,I+2) + u0(J+2,I-2) + u0(J+2,I+2)) / 4;\n",
+"\n",
+"for j = [J-1 J+1]\n",
+" for i = [I-1 I+1]\n",
+" u(j,i) = (u(j-1,i-1) + u(j-1,i+1) + u(j+1,i-1) + u(j+1,i+1)) / 4;\n",
+" end\n",
+"end\n",
+"\n",
+"j1 = [J-1 J J J+1];\n",
+"i1 = [I I-1 I+1 I];\n",
+"for k = 1:4\n",
+" i = i1(k);\n",
+" j = j1(k);\n",
+" u(j,i) = (u(j,i-1) + u(j,i+1) + u(j-1,i) + u(j+1,i)) / 4;\n",
+"end\n",
+"\n",
+"kf = 2;\n",
+"tab = zeros(kf+1,(m-2)*(n-2));\n",
+"row = [];\n",
+"for j = 2:n-1\n",
+" row = [row u(j,2:m-1)];\n",
+"end\n",
+"tab(1,:) = row;\n",
+"for k = 1:kf\n",
+" row = [];\n",
+" for j = 2:n-1\n",
+" for i = 2:m-1\n",
+" u(j,i) = (u(j,i-1) + u(j,i+1) + u(j-1,i) + u(j+1,i)) / 4;\n",
+" end\n",
+" row = [row u(j,2:m-1)];\n",
+" end\n",
+" row = round(row*10^4)/10^4;\n",
+" tab(k+1,:) = row;\n",
+"end\n",
+"mprintf('%4s %9s %9s %9s %9s %10s %10s %10s %10s %10s','r','u11','u21','u31','u12','u22','u32','u13','u23','u33')\n",
+"disp([(1:k+1)' tab])"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 10.4: Poisson_Equation_using_Liebmann_Iterative_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 10.4\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"h = 1/3;\n",
+"x = 0:h:1;\n",
+"y = 0:h:1;\n",
+"m = length(y);\n",
+"n = length(x);\n",
+"u = zeros(m,n);\n",
+"u(m,2:n-1) = 1;\n",
+"\n",
+"kf = 5;\n",
+"tab = zeros(kf,(m-2)*(n-2));\n",
+"for k = 1:kf\n",
+" row = [];\n",
+" for j = 2:n-1\n",
+" for i = 2:m-1\n",
+" constant = 10/9* (5 + 1/9*(i-1)^2 + 1/9*(j-1)^2);\n",
+" u(j,i) = (u(j,i-1) + u(j,i+1) + u(j-1,i) + u(j+1,i) + constant) / 4;\n",
+" end\n",
+" row = [row u(j,2:m-1)];\n",
+" end\n",
+" row = round(row*10^4)/10^4;\n",
+" tab(k,:) = row;\n",
+"end\n",
+"mprintf('%4s %9s %9s %9s %9s','r','u11','u21','u12','u22')\n",
+"disp([(1:k)' tab])"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 10.5: Laplace_Equation_using_Liebmann_Over_Relaxation_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 10.5\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 0:4;\n",
+"y = 0:4;\n",
+"m = length(y);\n",
+"n = length(x);\n",
+"u = zeros(m,n);\n",
+"u(m,:) = x.^3;\n",
+"u(:,n) = 16*y';\n",
+"u0 = u;\n",
+"\n",
+"I = ceil(m/2);\n",
+"J = ceil(n/2);\n",
+"\n",
+"u(J,I) = (u0(J-2,I-2) + u0(J-2,I+2) + u0(J+2,I-2) + u0(J+2,I+2)) / 4;\n",
+"\n",
+"for j = [J-1 J+1]\n",
+" for i = [I-1 I+1]\n",
+" u(j,i) = (u(j-1,i-1) + u(j-1,i+1) + u(j+1,i-1) + u(j+1,i+1)) / 4;\n",
+" end\n",
+"end\n",
+"\n",
+"j1 = [J-1 J J J+1];\n",
+"i1 = [I I-1 I+1 I];\n",
+"for k = 1:4\n",
+" i = i1(k);\n",
+" j = j1(k);\n",
+" u(j,i) = (u(j,i-1) + u(j,i+1) + u(j-1,i) + u(j+1,i)) / 4;\n",
+"end\n",
+"disp(u,'u:')\n",
+"\n",
+"p = m-1;\n",
+"q = n-1;\n",
+"c = cos(%pi/p) + cos(%pi/q);\n",
+"w = 4/(2+sqrt(4-c^2));\n",
+"w = round(w*10^3)/10^3;\n",
+"\n",
+"kf = 10;\n",
+"tab = zeros(kf+1,(m-2)*(n-2));\n",
+"row = [];\n",
+"for j = 2:n-1\n",
+" row = [row u(j,2:m-1)];\n",
+"end\n",
+"tab(1,:) = row;\n",
+"for k = 1:kf\n",
+" row = [];\n",
+" for j = 2:n-1\n",
+" for i = 2:m-1\n",
+" u(j,i) = (u(j,i-1) + u(j,i+1) + u(j-1,i) + u(j+1,i)) *w/4 + (1-w)*u(j,i);\n",
+" end\n",
+" row = [row u(j,2:m-1)];\n",
+" end\n",
+" row = round(row*10^4)/10^4;\n",
+" tab(k+1,:) = row;\n",
+"end\n",
+"mprintf('\n\n%8s %9s %10s %10s %9s %10s %10s %9s %9s','u11','u21','u31','u12','u22','u32','u13','u23','u33')\n",
+"disp(tab)"
+ ]
+ }
+],
+"metadata": {
+ "kernelspec": {
+ "display_name": "Scilab",
+ "language": "scilab",
+ "name": "scilab"
+ },
+ "language_info": {
+ "file_extension": ".sce",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "scilab",
+ "version": "0.7.1"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/11-Hyperbolic_Partial_Differential_Equations.ipynb b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/11-Hyperbolic_Partial_Differential_Equations.ipynb
new file mode 100644
index 0000000..25da4e2
--- /dev/null
+++ b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/11-Hyperbolic_Partial_Differential_Equations.ipynb
@@ -0,0 +1,137 @@
+{
+"cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Chapter 11: Hyperbolic Partial Differential Equations"
+ ]
+ },
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 11.1: Initial_Value_Problem_using_Wave_Equation.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 11.1\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"delx = 1/8;\n",
+"delt = 1/8;\n",
+"x = 0:delx:1;\n",
+"t = 0:delt:1;\n",
+"m = length(x);\n",
+"n = length(t);\n",
+"u =zeros(n,m);\n",
+"u(1,:) = sin(%pi*x);\n",
+"N = 1/delx;\n",
+"r = delt/delx;\n",
+"\n",
+"for j = 2:n\n",
+" for i = 2:m-1\n",
+" if j == 2 then\n",
+" u(j,i) = (2*(1-r^2)*u(j-1,i) + r^2*(u(j-1,i-1) + u(j-1,i+1))) /2;\n",
+" else\n",
+" u(j,i) = 2*(1-r^2)*u(j-1,i) + r^2*(u(j-1,i-1) + u(j-1,i+1)) - u(j-2,i);\n",
+" end\n",
+" end\n",
+"end\n",
+"u = round(u*10^4)/10^4;\n",
+"mprintf('\n\n%6s %9s %9s %8s %s %8s %10s %10s %9s %7s %s','t','x = 0.0','x = 0.125','x = 0.25','x = 0.375','x = 0.5','x=0.625','x = 0.75','x=0.875','x = 1.0','n');\n",
+"disp([(0:1/8:1)' u (0:n-1)']);\n",
+"mprintf('\n\n');\n",
+"t = [1/2 1];\n",
+"for i = 1:length(t)\n",
+" Ex(i,:) = sin(%pi*x) * cos(%pi*t(i));\n",
+"end\n",
+"Ex = round(Ex*10^4)/10^4;\n",
+"disp('At t = 1/2:')\n",
+"disp(u(find(x==1/2),:),'Computed Solution:')\n",
+"disp(Ex(1,:),'Actual Solution:')\n",
+"\n",
+"disp('At t = 1:')\n",
+"disp(u(find(x==1),:),'Computed Solution:')\n",
+"disp(Ex(2,:),'Actual Solution:')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 11.2: Initial_Value_Problem_using_Wave_Equation.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 11.2\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"delx = 0.2;\n",
+"delt = 0.2;\n",
+"x = 0:delx:1;\n",
+"t = 0:delt:0.8;\n",
+"m = length(x);\n",
+"n = length(t);\n",
+"u =zeros(n,m);\n",
+"u(1,:) = x^2;\n",
+"u(:,m) = 1+t';\n",
+"N = 1/delx;\n",
+"r = delt/delx;\n",
+"\n",
+"for j = 2:n\n",
+" for i = 2:m-1\n",
+" if j == 2 then\n",
+" u(j,i) = (2*(1-r^2)*u(j-1,i) + r^2*(u(j-1,i-1) + u(j-1,i+1)) + 2*delt) /2;\n",
+" else\n",
+" u(j,i) = 2*(1-r^2)*u(j-1,i) + r^2*(u(j-1,i-1) + u(j-1,i+1)) - u(j-2,i);\n",
+" end\n",
+" end\n",
+"end\n",
+"u = round(u*10^4)/10^4;\n",
+"mprintf('\n%5s %9s %7s %7s %s %6s %6s','t','x = 0.0','x = 0.2','x = 0.4','x = 0.6','x = 0.8','x = 1.0');\n",
+"disp([t' u])"
+ ]
+ }
+],
+"metadata": {
+ "kernelspec": {
+ "display_name": "Scilab",
+ "language": "scilab",
+ "name": "scilab"
+ },
+ "language_info": {
+ "file_extension": ".sce",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "scilab",
+ "version": "0.7.1"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/2-Solution_of_Algebraic_and_Transcendental_Equations.ipynb b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/2-Solution_of_Algebraic_and_Transcendental_Equations.ipynb
new file mode 100644
index 0000000..54c3da6
--- /dev/null
+++ b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/2-Solution_of_Algebraic_and_Transcendental_Equations.ipynb
@@ -0,0 +1,642 @@
+{
+"cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Chapter 2: Solution of Algebraic and Transcendental Equations"
+ ]
+ },
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 2.10: Newton_Formula.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 2.10\n",
+"clc\n",
+"clear\n",
+"\n",
+"N = 12;\n",
+"xold = 3.5;\n",
+"iter = 1;\n",
+"maxit = 3;\n",
+"\n",
+"while (1)\n",
+" xnew = (xold + N/xold) / 2;\n",
+" if iter == maxit then\n",
+" break\n",
+" end\n",
+" xold = xnew;\n",
+" iter = iter + 1;\n",
+"end\n",
+"root = round(xnew*10^4) / 10^4;\n",
+"disp(root,'root = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 2.11: Newton_Raphson_Extended_Formula.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"// Example 2.11\n",
+"// This is an analytical problem and need not be coded."
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 2.12: Root_using_Muller_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 2.12\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f] = fun12(x)\n",
+" f = x.^3 - x - 1;\n",
+"endfunction\n",
+"\n",
+"x = [0 1 2];\n",
+"h = [x(2)-x(1) x(3)-x(2)];\n",
+"lamdai = h(2)/h(1);\n",
+"deli = 1 + lamdai;\n",
+"f = fun12(x);\n",
+"\n",
+"g = f(1)*lamdai^2 - f(2)*deli^2 + f(3)*(lamdai + deli);\n",
+"lamda = -2*f(3)*deli / (g + sqrt(g^2 - 4*f(3)*deli*(f(1)*lamdai - f(2)*deli + f(3))));\n",
+"xnew = x(3) + lamda*h(2);\n",
+"xnew = round(xnew*10^5) / 10^5;\n",
+"disp(xnew,'root = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 2.13: Graeffe_Root_Squaring_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 2.13\n",
+"clc\n",
+"clear\n",
+"\n",
+"a = [-6 11 -6 1];\n",
+"maxit = 3;\n",
+"for iter = 1:maxit\n",
+" a = [a(4)^2 -(a(3)^2 -2*a(2)*a(4)) (a(2)^2 - 2*a(1)*a(3)) -a(1)^2];\n",
+" root = abs([a(4)/a(3) a(3)/a(2) a(2)/a(1)])^(1/(2^iter));\n",
+"end\n",
+"root = round(root*10^5) / 10^5;\n",
+"disp(root,'Estimated roots for the polynomial are: ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 2.1: Root_using_Bisection_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 2.1\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [root] = Bisection(fun,x,tol,maxit)\n",
+"// Bisection: Computes roots of the function in the given range using Bisection Method\n",
+"//// Input: Bisection(fun,x,tol,maxit)\n",
+"// fun = function handle\n",
+"// x = range in between sign change is evident\n",
+"// tol = Maximum error between iterations that can be tolerated\n",
+"// maxit = Maximum number of iterations\n",
+"//// Output: [root]\n",
+"// Root: Root of the given function in defined range\n",
+"\n",
+"if fun(x(1)) > 0 then\n",
+" xu = x(1); xl = x(2);\n",
+"else\n",
+" xu = x(2); xl = x(1);\n",
+"end\n",
+"\n",
+"Ea = 1;\n",
+"iter = 1;\n",
+"\n",
+"while(1)\n",
+" xr(iter) = (xl(iter) + xu(iter)) / 2;\n",
+" if fun(xr(iter)) > 0 then\n",
+" xu(iter+1) = xr(iter);\n",
+" xl(iter+1) = xl(iter);\n",
+" elseif fun(xr(iter)) < 0 then\n",
+" xl(iter+1) = xr(iter);\n",
+" xu(iter+1) = xu(iter);\n",
+" else\n",
+" break\n",
+" end\n",
+" \n",
+" if iter>1 then\n",
+" Ea(iter) = 100 * abs((xr(iter) - xr(iter-1)) / xr(iter));\n",
+" end\n",
+" \n",
+" if Ea(iter) < tol | iter == maxit then\n",
+" break\n",
+" end\n",
+" iter = iter + 1;\n",
+"end\n",
+"root = xr(iter);\n",
+"endfunction\n",
+"\n",
+"function f = fun1(x)\n",
+" f = x.^3 -9*x + 1;\n",
+"endfunction\n",
+"\n",
+"x = [2 4];\n",
+"tol = 1e-4;\n",
+"maxit = 5;\n",
+"root = Bisection(fun1,x,tol,maxit);\n",
+"disp(root,'root = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 2.2: Root_using_Regula_Falsi_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 2.2\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [root] = FalsePosition(fun,x,tol,maxit)\n",
+"// FalsePosition: Computes roots of the function in the given range using False Position Method\n",
+"//// Input: FalsePosition(fun,x,tol,maxit)\n",
+"// fun = function handle\n",
+"// x = range in between sign change is evident\n",
+"// tol = Maximum error between iterations that can be tolerated\n",
+"// maxit = Maximum number of iterations\n",
+"//// Output: [root]\n",
+"// Root: Root of the given function in defined range\n",
+"\n",
+"if fun(x(1)) > 0 then\n",
+" xu = x(1); xl = x(2);\n",
+"else\n",
+" xu = x(2); xl = x(1);\n",
+"end\n",
+"\n",
+"Ea = 1;\n",
+"iter = 1;\n",
+"\n",
+"while(1)\n",
+" xr(iter) = xl(iter) - ((xu(iter)-xl(iter)) / (fun(xu(iter))-fun(xl(iter))) * fun(xl(iter)));\n",
+" if fun(xr(iter)) > 0 then\n",
+" xu(iter+1) = xr(iter);\n",
+" xl(iter+1) = xl(iter);\n",
+" elseif fun(xr(iter)) < 0 then\n",
+" xl(iter+1) = xr(iter);\n",
+" xu(iter+1) = xu(iter);\n",
+" else\n",
+" break\n",
+" end\n",
+" \n",
+" if iter>1 then\n",
+" Ea(iter) = 100 * abs((xr(iter) - xr(iter-1)) / xr(iter));\n",
+" end\n",
+" \n",
+" if Ea(iter) < tol | iter == maxit then\n",
+" break\n",
+" end\n",
+" iter = iter + 1;\n",
+"end\n",
+"root = xr(iter);\n",
+"endfunction\n",
+"\n",
+"function f = fun1(x)\n",
+" f = x.^3 -9*x + 1;\n",
+"endfunction\n",
+"\n",
+"x = [2 4; 2 3];\n",
+"tol = 1e-4;\n",
+"maxit = 3;\n",
+"for i = 1:2\n",
+" root = FalsePosition(fun1,x(i,:),tol,maxit);\n",
+" root = round(root*10^5)/10^5;\n",
+" disp(strcat(['root(',string(i),') = ',string(root)]))\n",
+"end"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 2.3: Root_using_Regula_Falsi_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 2.3\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [root] = FalsePosition(fun,x,tol,maxit)\n",
+"// FalsePosition: Computes roots of the function in the given range using False Position Method\n",
+"//// Input: FalsePosition(fun,x,tol,maxit)\n",
+"// fun = function handle\n",
+"// x = range in between sign change is evident\n",
+"// tol = Maximum error between iterations that can be tolerated\n",
+"// maxit = Maximum number of iterations\n",
+"//// Output: [root]\n",
+"// Root: Root of the given function in defined range\n",
+"\n",
+"if fun(x(1)) > 0 then\n",
+" xu = x(1); xl = x(2);\n",
+"else\n",
+" xu = x(2); xl = x(1);\n",
+"end\n",
+"\n",
+"Ea = 1;\n",
+"iter = 1;\n",
+"\n",
+"while(1)\n",
+" xr(iter) = xl(iter) - ((xu(iter)-xl(iter)) / (fun(xu(iter))-fun(xl(iter))) * fun(xl(iter)));\n",
+" if fun(xr(iter)) > 0 then\n",
+" xu(iter+1) = xr(iter);\n",
+" xl(iter+1) = xl(iter);\n",
+" elseif fun(xr(iter)) < 0 then\n",
+" xl(iter+1) = xr(iter);\n",
+" xu(iter+1) = xu(iter);\n",
+" else\n",
+" break\n",
+" end\n",
+" \n",
+" if iter>1 then\n",
+" Ea(iter) = 100 * abs((xr(iter) - xr(iter-1)) / xr(iter));\n",
+" end\n",
+" \n",
+" if Ea(iter) < tol | iter == maxit then\n",
+" break\n",
+" end\n",
+" iter = iter + 1;\n",
+"end\n",
+"root = xr(iter);\n",
+"endfunction\n",
+"\n",
+"function f = fun3(x)\n",
+" f = log(x) - cos(x);\n",
+"endfunction\n",
+"\n",
+"x = [1 2];\n",
+"tol = 1e-4;\n",
+"maxit = 5;\n",
+"root = FalsePosition(fun3,x,tol,maxit);\n",
+"disp(round(root*10^4)/10^4,'root = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 2.4: Root_using_Regula_Falsi_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 2.4\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [root] = FalsePosition(fun,x,tol,maxit)\n",
+"// FalsePosition: Computes roots of the function in the given range using False Position Method\n",
+"//// Input: FalsePosition(fun,x,tol,maxit)\n",
+"// fun = function handle\n",
+"// x = range in between sign change is evident\n",
+"// tol = Maximum error between iterations that can be tolerated\n",
+"// maxit = Maximum number of iterations\n",
+"//// Output: [root]\n",
+"// Root: Root of the given function in defined range\n",
+"\n",
+"if fun(x(1)) > 0 then\n",
+" xu = x(1); xl = x(2);\n",
+"else\n",
+" xu = x(2); xl = x(1);\n",
+"end\n",
+"\n",
+"Ea = 1;\n",
+"iter = 1;\n",
+"\n",
+"while(1)\n",
+" xr(iter) = xl(iter) - ((xu(iter)-xl(iter)) / (fun(xu(iter))-fun(xl(iter))) * fun(xl(iter)));\n",
+" if fun(xr(iter)) > 0 then\n",
+" xu(iter+1) = xr(iter);\n",
+" xl(iter+1) = xl(iter);\n",
+" elseif fun(xr(iter)) < 0 then\n",
+" xl(iter+1) = xr(iter);\n",
+" xu(iter+1) = xu(iter);\n",
+" else\n",
+" break\n",
+" end\n",
+" \n",
+" if iter>1 then\n",
+" Ea(iter) = 100 * abs((xr(iter) - xr(iter-1)) / xr(iter));\n",
+" end\n",
+" \n",
+" if Ea(iter) < tol | iter == maxit then\n",
+" break\n",
+" end\n",
+" iter = iter + 1;\n",
+"end\n",
+"root = xr(iter);\n",
+"endfunction\n",
+"\n",
+"function f = fun4(x)\n",
+" f = x.*log10(x) - 1.2;\n",
+"endfunction\n",
+"\n",
+"clc\n",
+"x = [2 3];\n",
+"tol = 1e-4;\n",
+"maxit = 2;\n",
+"root = FalsePosition(fun4,x,tol,maxit);\n",
+"disp(round(root*10^4)/10^4,'root = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 2.5: Root_using_Method_of_Iteration.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 2.5\n",
+"clc\n",
+"clear\n",
+"\n",
+"function f = fun5(x)\n",
+" f = exp(-x)/10;\n",
+"endfunction\n",
+"\n",
+"clc\n",
+"tol = 1e-4;\n",
+"maxit = 4;\n",
+"xold = 0;\n",
+"iter = 1;\n",
+"while(1)\n",
+" xnew = fun5(xold);\n",
+" EA = abs(xnew - xold);\n",
+" if EA < tol | iter > maxit then\n",
+" break\n",
+" end\n",
+" xold = xnew;\n",
+" iter = iter + 1;\n",
+"end\n",
+"root = round(xnew*10^4) / 10^4; //rounded to 4 decimal places\n",
+"disp(root,'root = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 2.6: Root_using_Method_of_Iteration.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 2.6\n",
+"clc\n",
+"clear\n",
+"\n",
+"function f = fun6(x)\n",
+" f = 1./ sqrt(x+1);\n",
+"endfunction\n",
+"\n",
+"tol = 1e-4;\n",
+"maxit = 6;\n",
+"xold = 1;\n",
+"iter = 1;\n",
+"while(1)\n",
+" xnew = fun6(xold);\n",
+" EA = abs(xnew - xold);\n",
+" if EA < tol | iter > maxit then\n",
+" break\n",
+" end\n",
+" xold = xnew;\n",
+" iter = iter + 1;\n",
+"end\n",
+"root = round(xnew*10^4) / 10^4; //rounded to 4 decimal places\n",
+"disp(root,'root = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 2.7: Root_using_Newton_Raphson_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 2.7\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f,df] = fun7(x)\n",
+" f = x.*exp(x) - 2;\n",
+" df = x.*exp(x) + exp(x);\n",
+"endfunction\n",
+"\n",
+"xold = 1;\n",
+"maxit = 2;\n",
+"iter = 1;\n",
+"\n",
+"while (1)\n",
+" [fx,dfx] = fun7(xold);\n",
+" xnew = xold - fx/dfx;\n",
+" if iter == maxit then\n",
+" break\n",
+" end\n",
+" xold = xnew;\n",
+" iter = iter + 1;\n",
+"end\n",
+"root = round(xnew*10^3) / 10^3;\n",
+"disp(root,'root = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 2.8: Root_using_Newton_Raphson_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 2.8\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f,df] = fun8(x)\n",
+" f = x.^3 - x - 1;\n",
+" df = 3*x.^2 - 1;\n",
+"endfunction\n",
+"\n",
+"xold = 1;\n",
+"maxit = 5;\n",
+"iter = 1;\n",
+"\n",
+"while (1)\n",
+" [fx,dfx] = fun8(xold);\n",
+" xnew = xold - fx/dfx;\n",
+" if iter == maxit then\n",
+" break\n",
+" end\n",
+" xold = xnew;\n",
+" iter = iter + 1;\n",
+"end\n",
+"root = round(xnew*10^4) / 10^4;\n",
+"disp(root,'root = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 2.9: Newton_Scheme_of_Iteration.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"// Example 2.9\n",
+"// This is an analytical problem and need not be coded."
+ ]
+ }
+],
+"metadata": {
+ "kernelspec": {
+ "display_name": "Scilab",
+ "language": "scilab",
+ "name": "scilab"
+ },
+ "language_info": {
+ "file_extension": ".sce",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "scilab",
+ "version": "0.7.1"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/3-Solution_of_Linear_System_of_Equations_and_Matrix_Inversion.ipynb b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/3-Solution_of_Linear_System_of_Equations_and_Matrix_Inversion.ipynb
new file mode 100644
index 0000000..4450d1a
--- /dev/null
+++ b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/3-Solution_of_Linear_System_of_Equations_and_Matrix_Inversion.ipynb
@@ -0,0 +1,549 @@
+{
+"cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Chapter 3: Solution of Linear System of Equations and Matrix Inversion"
+ ]
+ },
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 3.10: Matrix_Inverse_using_Gauss_Jordan_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 3.10\n",
+"clc\n",
+"clear\n",
+"\n",
+"A = [1 1 1; 4 3 -1; 3 5 3];\n",
+"n = length (A(1,:));\n",
+"Aug = [A,eye(n,n)];\n",
+"\n",
+"N = 1:n;\n",
+"for i = 1:n\n",
+" dummy1 = N;\n",
+" dummy1(i) = [];\n",
+" index(i,:) = dummy1;\n",
+"end\n",
+"\n",
+"// Forward Elimination\n",
+"for j = 1:n\n",
+" [dummy2,t] = max(abs(Aug(j:n,j)));\n",
+" lrow = t+j-1;\n",
+" Aug([j,lrow],:) = Aug([lrow,j],:);\n",
+" Aug(j,:) = Aug(j,:) / Aug(j,j);\n",
+" for i = index(j,:)\n",
+" Aug(i,:) = Aug(i,:) - Aug(i,j) / Aug(j,j) * Aug(j,:);\n",
+" end\n",
+"end\n",
+"Inv_A = Aug(:,n+1:2*n);\n",
+"disp(Inv_A,'Inverse of A (A-1) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 3.1: Gauss_Elimination_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 3.1\n",
+"clc\n",
+"clear\n",
+"\n",
+"A = [2 3 -1; 4 4 -3; -2 3 -1]; //Coefficient Matrix\n",
+"B = [5; 3; 1]; //Constant Matrix\n",
+"\n",
+"n = length(B);\n",
+"Aug = [A,B];\n",
+"\n",
+"// Forward Elimination\n",
+"for j = 1:n-1\n",
+" for i = j+1:n\n",
+" Aug(i,j:n+1) = Aug(i,j:n+1) - Aug(i,j) / Aug(j,j) * Aug(j,j:n+1);\n",
+" end\n",
+"end\n",
+"\n",
+"// Backward Substitution\n",
+"x = zeros(n,1);\n",
+"x(n) = Aug(n,n+1) / Aug(n,n);\n",
+"for i = n-1:-1:1\n",
+" x(i) = (Aug(i,n+1)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);\n",
+"end\n",
+"disp(strcat(['x = ',string(x(1))]))\n",
+"disp(strcat(['y = ',string(x(2))]))\n",
+"disp(strcat(['z = ',string(x(3))]))\n",
+""
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 3.2: Gauss_Elimination_Method_with_Partial_Pivoting.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 3.2\n",
+"clc\n",
+"clear\n",
+"\n",
+"A = [1 1 1; 3 3 4; 2 1 3]; //Coefficient Matrix\n",
+"B = [7; 24; 16]; //Constant Matrix\n",
+"\n",
+"n = length(B);\n",
+"Aug = [A,B];\n",
+"\n",
+"// Forward Elimination\n",
+"for j = 1:n-1\n",
+" // Partial Pivoting\n",
+" [dummy,t] = max(abs(Aug(j:n,j)));\n",
+" lrow = t(1)+j-1;\n",
+" Aug([j,lrow],:) = Aug([lrow,j],:);\n",
+" \n",
+" for i = j+1:n\n",
+" Aug(i,j:n+1) = Aug(i,j:n+1) - Aug(i,j) / Aug(j,j) * Aug(j,j:n+1);\n",
+" end\n",
+"end\n",
+"\n",
+"// Backward Substitution\n",
+"x = zeros(n,1);\n",
+"x(n) = Aug(n,n+1) / Aug(n,n);\n",
+"for i = n-1:-1:1\n",
+" x(i) = (Aug(i,n+1)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);\n",
+"end\n",
+"disp(strcat(['x = ',string(x(1))]))\n",
+"disp(strcat(['y = ',string(x(2))]))\n",
+"disp(strcat(['z = ',string(x(3))]))\n",
+""
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 3.3: Gauss_Elimination_Method_with_Partial_Pivoting.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 3.3\n",
+"clc\n",
+"clear\n",
+"\n",
+"A = [0 4 2 8; 4 10 5 4; 4 5 6.5 2; 9 4 4 0]; //Coefficient Matrix\n",
+"B = [24; 32; 26; 21]; //Constant Matrix\n",
+"\n",
+"n = length(B);\n",
+"Aug = [A,B];\n",
+"\n",
+"// Forward Elimination\n",
+"for j = 1:n-1\n",
+" // Partial Pivoting\n",
+" [dummy,t] = max(abs(Aug(j:n,j)));\n",
+" lrow = t(1)+j-1;\n",
+" Aug([j,lrow],:) = Aug([lrow,j],:);\n",
+" \n",
+" for i = j+1:n\n",
+" Aug(i,j:n+1) = Aug(i,j:n+1) - Aug(i,j) / Aug(j,j) * Aug(j,j:n+1);\n",
+" end\n",
+"end\n",
+"\n",
+"// Backward Substitution\n",
+"x = zeros(n,1);\n",
+"x(n) = Aug(n,n+1) / Aug(n,n);\n",
+"for i = n-1:-1:1\n",
+" x(i) = (Aug(i,n+1)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);\n",
+"end\n",
+"\n",
+"disp(strcat(['x1 = ',string(x(1))]))\n",
+"disp(strcat(['x2 = ',string(x(2))]))\n",
+"disp(strcat(['x3 = ',string(x(3))]))\n",
+"disp(strcat(['x4 = ',string(x(4))]))"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 3.4: Gauss_Jordan_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 3.4\n",
+"clc\n",
+"clear\n",
+"\n",
+"A = [1 2 1; 2 3 4; 4 3 2];\n",
+"B = [8; 20; 16];\n",
+"n = length (B);\n",
+"Aug = [A,B];\n",
+"\n",
+"// Forward Elimination\n",
+"for j = 1:n-1\n",
+" for i = j+1:n\n",
+" Aug(i,j:n+1) = Aug(i,j:n+1) - Aug(i,j) / Aug(j,j) * Aug(j,j:n+1);\n",
+" end\n",
+"end\n",
+"\n",
+"// Backward Elimination\n",
+"for j = n:-1:2\n",
+" Aug(1:j-1,:) = Aug(1:j-1,:) - Aug(1:j-1,j) / Aug(j,j) * Aug(j,:);\n",
+"end\n",
+"\n",
+"// Diagonal Normalization\n",
+"for j=1:n\n",
+" Aug(j,:) = Aug(j,:) / Aug(j,j);\n",
+"end\n",
+"x = Aug(:,n+1);\n",
+"disp(strcat(['x = ',string(x(1))]))\n",
+"disp(strcat(['y = ',string(x(2))]))\n",
+"disp(strcat(['z = ',string(x(3))]))"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 3.5: Crout_Reduction_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 3.5\n",
+"clc\n",
+"clear\n",
+"\n",
+"A = [5 -2 1; 7 1 -5; 3 7 4];\n",
+"B = [4; 8; 10];\n",
+"\n",
+"n = length (B);\n",
+"L = zeros(n,n); // L = Lower Triangular Matrix Initiation\n",
+"U = eye(n,n); // U = Upper Triangular Matrix Initiation\n",
+"\n",
+"// LU Decomposition\n",
+"for i = 1:n\n",
+" sum1 = zeros(n-i+1,1);\n",
+" for k = 1:i-1\n",
+" sum1 = sum1 + L(i:n,k) * U(k,i);\n",
+" end\n",
+" L(i:n,i) = A(i:n,i) - sum1;\n",
+"\n",
+" sum2 = zeros(1,n-i);\n",
+" for k = 1:i-1\n",
+" sum2 = sum2 + L(i,k) * U(k,i+1:n);\n",
+" end\n",
+" U(i,i+1:n) = (A(i,i+1:n) - sum2) / L(i,i);\n",
+"end\n",
+"\n",
+"// Forward Substitution\n",
+"D = ones(n,1);\n",
+"for i = 1:n\n",
+" sum3 = 0;\n",
+" for k = 1:i-1\n",
+" sum3 = sum3 + L(i,k) * D(k);\n",
+" end\n",
+" D(i) = (B(i) - sum3) / L(i,i);\n",
+"end\n",
+"\n",
+"// Back Substitution\n",
+"x = ones(n,1);\n",
+"for i = n:-1:1\n",
+" sum4 = 0;\n",
+" for k = i+1:n\n",
+" sum4 = sum4 + U(i,k) * x(k);\n",
+" end\n",
+" x(i) = D(i) - sum4;\n",
+"end\n",
+"\n",
+"disp(strcat(['x1 = ',string(x(1))]))\n",
+"disp(strcat(['x2 = ',string(x(2))]))\n",
+"disp(strcat(['x3 = ',string(x(3))]))\n",
+""
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 3.6: Jacobi_Iterative_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 3.6\n",
+"clc\n",
+"clear\n",
+"\n",
+"A = [83 11 -4; 7 52 13; 3 8 29];\n",
+"B = [95; 104; 71];\n",
+"\n",
+"n = length (B);\n",
+"tol = 1e-4;\n",
+"iter = 1;\n",
+"maxit = 5;\n",
+"\n",
+"x = zeros(n,1); //Intial guess\n",
+"E = ones(n,1); //Assuming to avoid variable size error\n",
+"S = diag(diag(A));\n",
+"\n",
+"\n",
+"while (1)\n",
+" x(:,iter+1) = S\(B + (S-A)*(x(:,iter)));\n",
+" E(:,iter+1) = (x(:,iter+1)-x(:,iter))./x(:,iter+1)*100;\n",
+" if x(:,iter) == 0\n",
+" Error = 1;\n",
+" else\n",
+" Error = sqrt((sum((E(:,iter+1)).^2))/n);\n",
+" end\n",
+" \n",
+" if Error <= tol | iter == maxit\n",
+" break\n",
+" end\n",
+" iter = iter+1;\n",
+"end\n",
+"xact = x(:,iter);\n",
+"x = round(x*10^4)/10^4;\n",
+"x(:,1) = [];\n",
+"mprintf('%s %3s %9s %9s','Iter No.','x','y','z');\n",
+"disp([(1:iter)' x']);"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 3.7: Gauss_Seidel_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 3.7\n",
+"clc\n",
+"clear\n",
+"\n",
+"A = [1 -1/4 -1/4 0; -1/4 1 0 -1/4; -1/4 0 1 -1/4; 0 -1/4 -1/4 1];\n",
+"B = [1/2; 1/2; 1/4; 1/4];\n",
+"\n",
+"n = length (B);\n",
+"tol = 1e-4;\n",
+"iter = 1;\n",
+"maxit = 5;\n",
+"\n",
+"x = zeros(n,1); //Intial guess\n",
+"E = ones(n,1); //Assuming to avoid variable size error\n",
+"S = diag(diag(A));\n",
+"T = S-A;\n",
+"xold = x;\n",
+"\n",
+"while (1)\n",
+" for i = 1:n\n",
+" x(i,iter+1) = (B(i) + T(i,:) * xold) / A(i,i);\n",
+" E(i,iter+1) = (x(i,iter+1)-xold(i))/x(i,iter+1)*100;\n",
+" xold(i) = x(i,iter+1);\n",
+" end\n",
+" \n",
+" if x(:,iter) == 0\n",
+" E = 1;\n",
+" else\n",
+" E = sqrt((sum((E(:,iter+1)).^2))/n);\n",
+" end\n",
+" \n",
+" if E <= tol | iter == maxit\n",
+" break\n",
+" end\n",
+" iter = iter + 1;\n",
+"end\n",
+"X = x(:,iter);\n",
+"x = round(x*10^5)/10^5;\n",
+"x(:,1) = [];\n",
+"mprintf('%s %3s %11s %10s %10s','Iter No.','x1','x2','x3','x4');\n",
+"disp([(1:iter)' x']);"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 3.8: Relaxation_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 3.8\n",
+"clc\n",
+"clear\n",
+"\n",
+"A = [6 -3 1; 2 1 -8;1 -7 1];\n",
+"b = [11; -15; 10];\n",
+"\n",
+"n = length (b);\n",
+"tol = 1e-4;\n",
+"iter = 1;\n",
+"maxit = 9;\n",
+"\n",
+"x = zeros(1,n); //Intial guess\n",
+"absA = abs(A);\n",
+"[dummy,index] = max(absA(1,:),absA(2,:),absA(3,:));\n",
+"if length(unique(index)) == n\n",
+" nu_T = diag(diag(A(index,:))) - A(index,:);\n",
+" if abs(diag(A(index,:))) - (sum(abs(nu_T),2)) > 0\n",
+" A(index,:) = A;\n",
+" b(index,:) = b;\n",
+" end\n",
+"end\n",
+"\n",
+"for iter = 1:maxit\n",
+" R(iter,:) = b' - x(iter,:) * A';\n",
+" [mx,i] = max(abs(R(iter,:)));\n",
+" Rmax(iter) = R(iter,i);\n",
+" dx(iter) = Rmax(iter) ./ A(i,i);\n",
+" x(iter+1,:) = x(iter,:);\n",
+" x(iter+1,i) = x(iter,i) + dx(iter);\n",
+"end\n",
+"R = round(R*10^4)/10^4;\n",
+"Rmax = round(Rmax*10^4)/10^4;\n",
+"dx = round(dx*10^4)/10^4;\n",
+"x = round(x*10^4)/10^4;\n",
+"mprintf('%s %3s %9s %9s %12s %10s %6s %9s %9s','Iter No.','R1','R2','R3','Max Ri','Diff dxi','x1','x2','x3');\n",
+"disp([(1:maxit)' R Rmax dx x(1:maxit,:)])"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 3.9: Matrix_Inverse_using_Gauss_Elimination_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 3.9\n",
+"clc\n",
+"clear\n",
+"\n",
+"A = [1 1 1; 4 3 -1; 3 5 3];\n",
+"n = length (A(1,:));\n",
+"Aug = [A,eye(n,n)];\n",
+"\n",
+"// Forward Elimination\n",
+"for j = 1:n-1\n",
+" for i = j+1:n\n",
+" Aug(i,j:2*n) = Aug(i,j:2*n) - Aug(i,j) / Aug(j,j) * Aug(j,j:2*n);\n",
+" end\n",
+"end\n",
+"\n",
+"// Backward Elimination\n",
+"for j = n:-1:2\n",
+" Aug(1:j-1,:) = Aug(1:j-1,:) - Aug(1:j-1,j) / Aug(j,j) * Aug(j,:);\n",
+"end\n",
+"\n",
+"// Diagonal Normalization\n",
+"for j=1:n\n",
+" Aug(j,:) = Aug(j,:) / Aug(j,j);\n",
+"end\n",
+"Inv_A = Aug(:,n+1:2*n);\n",
+"disp(Inv_A,'Inverse of A (A-1) = ')"
+ ]
+ }
+],
+"metadata": {
+ "kernelspec": {
+ "display_name": "Scilab",
+ "language": "scilab",
+ "name": "scilab"
+ },
+ "language_info": {
+ "file_extension": ".sce",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "scilab",
+ "version": "0.7.1"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/4-Eigenvalue_Problems.ipynb b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/4-Eigenvalue_Problems.ipynb
new file mode 100644
index 0000000..89624e0
--- /dev/null
+++ b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/4-Eigenvalue_Problems.ipynb
@@ -0,0 +1,278 @@
+{
+"cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Chapter 4: Eigenvalue Problems"
+ ]
+ },
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 4.1: Eigenvalues_and_Eigenvectors.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 4.1\n",
+"clc\n",
+"clear\n",
+"\n",
+"A = [2 3 2; 4 3 5; 3 2 9];\n",
+"v = [1; 1; 1];\n",
+"iter = 1;\n",
+"maxit = 5;\n",
+"\n",
+"while(1)\n",
+" u(:,iter) = A * v(:,iter);\n",
+" q(iter) = max(u(:,iter));\n",
+" v(:,iter+1) = u(:,iter) / q(iter);\n",
+" if iter == maxit then\n",
+" break\n",
+" end\n",
+" iter = iter + 1;\n",
+"end\n",
+"X = round(v(:,iter)*10^2) / 10^2;\n",
+"disp(X,'Eigen Vector:')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 4.2: Eigenvalues_and_Eigenvectors_using_Jacobi_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 4.2\n",
+"clc\n",
+"clear\n",
+"\n",
+"rt2 = sqrt(2);\n",
+"A = [1 rt2 2; rt2 3 rt2; 2 rt2 1];\n",
+"[n,n] = size(A);\n",
+"iter = 1;\n",
+"maxit = 3;\n",
+"D = A;\n",
+"S = 1;\n",
+"\n",
+"while(1)\n",
+" D_offdiag = D - diag(diag(D));\n",
+" [mx,index1] = max(abs(D_offdiag));\n",
+" i = index1(1);\n",
+" j = index1(2);\n",
+" if (D(i,i)-D(j,j)) == 0 then\n",
+" theta = %pi/4;\n",
+" else\n",
+" theta = atan(2*D(i,j)/(D(i,i)-D(j,j))) / 2;\n",
+" end\n",
+" S1 = eye(n,n);\n",
+" S1(i,i) = cos(theta);\n",
+" S1(i,j) = -sin(theta);\n",
+" S1(j,i) = sin(theta);\n",
+" S1(j,j) = cos(theta);\n",
+"\n",
+" D1 = inv(S1) * D * S1;\n",
+" for j = 1:n\n",
+" for i = 1:n\n",
+" if abs(D1(i,j)) < 1D-10 then\n",
+" D1(i,j) = 0;\n",
+" end\n",
+" end\n",
+" end\n",
+" S = S * S1;\n",
+"\n",
+" if D1 - diag(diag(D1)) == zeros(n,n) | iter == maxit then\n",
+" eigval = diag(D1);\n",
+" disp('Eigen Values:')\n",
+" disp(eigval)\n",
+"\n",
+" disp('Eigen Vectors:')\n",
+" disp(S(:,1))\n",
+" disp(S(:,2))\n",
+" disp(S(:,3))\n",
+" break\n",
+" end\n",
+"\n",
+" iter = iter + 1;\n",
+" D = D1;\n",
+"end"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 4.3: Eigenvalues_using_Jacobi_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 4.3\n",
+"clc\n",
+"clear\n",
+"\n",
+"A = [2 -1 0; -1 2 -1; 0 -1 2];\n",
+"[n,n] = size(A);\n",
+"iter = 1;\n",
+"maxit = 3;\n",
+"//Note: Diagonal form may be achieved at iter = 9. Modify maxit to greater than 9 for exact result.\n",
+"\n",
+"D = A;\n",
+"S = 1;\n",
+"\n",
+"while(1)\n",
+" D_offdiag = D - diag(diag(D));\n",
+" [mx,index1] = max(abs(D_offdiag));\n",
+" i = index1(1);\n",
+" j = index1(2);\n",
+" if (D(i,i)-D(j,j)) == 0 then\n",
+" theta = %pi/4;\n",
+" else\n",
+" theta = atan(2*D(i,j)/(D(i,i)-D(j,j))) / 2;\n",
+" end\n",
+" S1 = eye(n,n);\n",
+" S1(i,i) = cos(theta);\n",
+" S1(i,j) = -sin(theta);\n",
+" S1(j,i) = sin(theta);\n",
+" S1(j,j) = cos(theta);\n",
+"\n",
+" D1 = inv(S1) * D * S1;\n",
+" for j = 1:n\n",
+" for i = 1:n\n",
+" if abs(D1(i,j)) < 1D-10 then\n",
+" D1(i,j) = 0;\n",
+" end\n",
+" end\n",
+" end\n",
+" S = S * S1;\n",
+"\n",
+" if D1 - diag(diag(D1)) == zeros(n,n) | iter == maxit then\n",
+" eigval = diag(D1);\n",
+" eigval = round(eigval*10^3)/10^3;\n",
+" disp('Eigen Values:')\n",
+" disp(eigval)\n",
+" break\n",
+" end\n",
+"\n",
+" iter = iter + 1;\n",
+" D = D1;\n",
+"end"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 4.4: Eigenvalues_and_Eigenvectors_using_Jacobi_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 4.4\n",
+"clc\n",
+"clear\n",
+"\n",
+"A = [5 0 1; 0 -2 0; 1 0 5];\n",
+"[n,n] = size(A);\n",
+"iter = 1;\n",
+"maxit = 3;\n",
+"D = A;\n",
+"S = 1;\n",
+"\n",
+"while(1)\n",
+" D_offdiag = D - diag(diag(D));\n",
+" [mx,index1] = max(abs(D_offdiag));\n",
+" i = index1(1);\n",
+" j = index1(2);\n",
+" if (D(i,i)-D(j,j)) == 0 then\n",
+" theta = %pi/4;\n",
+" else\n",
+" theta = atan(2*D(i,j)/(D(i,i)-D(j,j))) / 2;\n",
+" end\n",
+" S1 = eye(n,n);\n",
+" S1(i,i) = cos(theta);\n",
+" S1(i,j) = -sin(theta);\n",
+" S1(j,i) = sin(theta);\n",
+" S1(j,j) = cos(theta);\n",
+"\n",
+" D1 = inv(S1) * D * S1;\n",
+" for j = 1:n\n",
+" for i = 1:n\n",
+" if abs(D1(i,j)) < 1D-10 then\n",
+" D1(i,j) = 0;\n",
+" end\n",
+" end\n",
+" end\n",
+" S = S * S1;\n",
+"\n",
+" if D1 - diag(diag(D1)) == zeros(n,n) | iter == maxit then\n",
+" eigval = diag(D1);\n",
+" disp('Eigen Values:')\n",
+" disp(eigval)\n",
+"\n",
+" disp('Eigen Vectors:')\n",
+" disp(S(:,1))\n",
+" disp(S(:,2))\n",
+" disp(S(:,3))\n",
+" break\n",
+" end\n",
+"\n",
+" iter = iter + 1;\n",
+" D = D1;\n",
+"end"
+ ]
+ }
+],
+"metadata": {
+ "kernelspec": {
+ "display_name": "Scilab",
+ "language": "scilab",
+ "name": "scilab"
+ },
+ "language_info": {
+ "file_extension": ".sce",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "scilab",
+ "version": "0.7.1"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/5-Curve_Fitting.ipynb b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/5-Curve_Fitting.ipynb
new file mode 100644
index 0000000..3be9134
--- /dev/null
+++ b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/5-Curve_Fitting.ipynb
@@ -0,0 +1,441 @@
+{
+"cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Chapter 5: Curve Fitting"
+ ]
+ },
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 5.10: Method_of_Moments.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 5.10\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 3:7;\n",
+"y = [31.9 34.6 33.8 27 31.6];\n",
+"\n",
+"delx = x(2) - x(1);\n",
+"mu1 = delx * sum(y);\n",
+"mu2 = delx * sum(x.*y);\n",
+"mu3 = delx * sum(x^2 .*y);\n",
+"\n",
+"n = length(y);\n",
+"l = x(1) - delx/2;\n",
+"u = x(n) + delx/2;\n",
+"\n",
+"t0 = u-l;\n",
+"t1 = integrate('x','x',l,u);\n",
+"t2 = integrate('x^2','x',l,u);\n",
+"t3 = integrate('x^3','x',l,u);\n",
+"t4 = integrate('x^4','x',l,u);\n",
+"\n",
+"M1 = [t2 t1 t0; t3 t2 t1; t4 t3 t2];\n",
+"M2 = [mu1; mu2; mu3];\n",
+"M1 = round(M1*10^2)/10^2;\n",
+"M = M1\M2;\n",
+"\n",
+"c = M(1);\n",
+"b = M(2);\n",
+"a = M(3);\n",
+"\n",
+"disp(round(a*10^4)/10^4, 'a =')\n",
+"disp(round(b*10^4)/10^4, 'b =')\n",
+"disp(round(c*10^4)/10^4, 'c =')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 5.1: Method_of_Group_Averages.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 5.1\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 10:10:80;\n",
+"y = [1.06 1.33 1.52 1.68 1.81 1.91 2.01 2.11];\n",
+"\n",
+"X = log(x);\n",
+"Y = log(y);\n",
+"\n",
+"n = length(Y);\n",
+"M1 = [sum(Y); sum(X.*Y)];\n",
+"M2 = [n sum(X); sum(X) sum(X.^2)];\n",
+"\n",
+"A = M2\M1;\n",
+"\n",
+"m = exp(A(1));\n",
+"n = A(2);\n",
+"\n",
+"disp(round(m*10^4)/10^4, 'm =')\n",
+"disp(round(n*10^4)/10^4, 'n =')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 5.2: Method_of_Group_Averages.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 5.2\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = [20 30 35 40 45 50];\n",
+"y = [10 11 11.8 12.4 13.5 14.4];\n",
+"\n",
+"X = x.^2;\n",
+"Y = y;\n",
+"\n",
+"n = length(Y);\n",
+"M1 = [sum(Y); sum(X.*Y)];\n",
+"M2 = [n sum(X); sum(X) sum(X.^2)];\n",
+"\n",
+"A = M2\M1;\n",
+"\n",
+"a = A(1);\n",
+"b = A(2);\n",
+"\n",
+"disp(round(a*10^4)/10^4, 'a =')\n",
+"disp(round(b*10^4)/10^4, 'b =')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 5.3: Method_of_Group_Averages.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 5.3\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = [8 10 15 20 30 40];\n",
+"y = [13 14 15.4 16.3 17.2 17.8];\n",
+"\n",
+"X = 1 ./x;\n",
+"Y = 1 ./y;\n",
+"\n",
+"n = length(Y);\n",
+"M1 = [sum(Y); sum(X.*Y)];\n",
+"M2 = [n sum(X); sum(X) sum(X.^2)];\n",
+"\n",
+"A = M2\M1;\n",
+"\n",
+"b = A(1);\n",
+"a = A(2);\n",
+"\n",
+"disp(round(a*10^4)/10^4, 'a =')\n",
+"disp(round(b*10^4)/10^4, 'b =')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 5.4: Method_of_Least_Squares.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 5.4\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"X = 0.5:0.5:3;\n",
+"Y = [15 17 19 14 10 7];\n",
+"\n",
+"n = length(Y);\n",
+"M1 = [sum(Y); sum(X.*Y)];\n",
+"M2 = [n sum(X); sum(X) sum(X.^2)];\n",
+"\n",
+"A = M2\M1;\n",
+"\n",
+"b = A(1);\n",
+"a = A(2);\n",
+"\n",
+"disp(round(a*10^4)/10^4, 'a =')\n",
+"disp(round(b*10^4)/10^4, 'b =')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 5.5: Method_of_Least_Squares.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 5.5\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 1:6;\n",
+"y = [2.6 5.4 8.7 12.1 16 20.2];\n",
+"\n",
+"X = x;\n",
+"Y = y ./x;\n",
+"\n",
+"n = length(Y);\n",
+"M1 = [sum(Y); sum(X.*Y)];\n",
+"M2 = [n sum(X); sum(X) sum(X.^2)];\n",
+"\n",
+"A = M2\M1;\n",
+"\n",
+"a = A(1);\n",
+"b = A(2);\n",
+"\n",
+"disp(round(a*10^5)/10^5, 'a =')\n",
+"disp(round(b*10^5)/10^5, 'b =')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 5.6: Method_of_Least_Squares.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 5.6\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"X = 1:0.2:2;\n",
+"Y = [0.98 1.4 1.86 2.55 2.28 3.2];\n",
+"\n",
+"n = length(Y);\n",
+"M1 = [sum(X.^4) sum(X.^3) sum(X.^2); sum(X.^3) sum(X.^2) sum(X); sum(X.^2) sum(X) n];\n",
+"M2 = [sum(X.^2 .* Y); sum(X.*Y); sum(Y)];\n",
+"A = M1\M2;\n",
+"\n",
+"a = A(1);\n",
+"b = A(2);\n",
+"c = A(3);\n",
+"\n",
+"disp(round(a*10^4)/10^4, 'a =')\n",
+"disp(round(b*10^4)/10^4, 'b =')\n",
+"disp(round(c*10^4)/10^4, 'c =')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 5.7: Method_of_Least_Squares.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 5.7\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 2:5;\n",
+"y = [27.8 62.1 110 161];\n",
+"\n",
+"X = log(x);\n",
+"Y = log(y);\n",
+"\n",
+"n = length(Y);\n",
+"M1 = [sum(X.^2) sum(X); sum(X) n];\n",
+"M2 = [sum(X.*Y); sum(Y)];\n",
+"M = M1\M2;\n",
+"\n",
+"b = M(1);\n",
+"A = M(2);\n",
+"a = exp(A);\n",
+"\n",
+"disp(round(a*10^4)/10^4, 'a =')\n",
+"disp(round(b*10^4)/10^4, 'b =')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 5.8: Principle_of_Least_Squares.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 5.8\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 1:4;\n",
+"y = [1.65 2.7 4.5 7.35];\n",
+"\n",
+"X = x;\n",
+"Y = log10(y);\n",
+"\n",
+"n = length(Y);\n",
+"M1 = [sum(X.^2) sum(X); sum(X) n];\n",
+"M2 = [sum(X.*Y); sum(Y)];\n",
+"M = M1\M2;\n",
+"\n",
+"B = M(1);\n",
+"A = M(2);\n",
+"a = 10^A;\n",
+"b = B/log10(%e);\n",
+"\n",
+"disp(round(a), 'a =')\n",
+"disp(round(b*10^4)/10^4, 'b =')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 5.9: Method_of_Moments.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 5.9\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 2:5;\n",
+"y = [27 40 55 68];\n",
+"\n",
+"delx = x(2) - x(1);\n",
+"mu1 = delx * sum(y);\n",
+"mu2 = delx * sum(x.*y);\n",
+"\n",
+"n = length(y);\n",
+"l = x(1) - delx/2;\n",
+"u = x(n) + delx/2;\n",
+"\n",
+"M1 = [integrate('x','x',l,u) u-l; integrate('x^2','x',l,u) integrate('x','x',l,u)];\n",
+"M2 = [mu1; mu2];\n",
+"M = M1\M2;\n",
+"\n",
+"a = M(1);\n",
+"b = M(2);\n",
+"\n",
+"disp(round(a*10^4)/10^4, 'a =')\n",
+"disp(round(b*10^4)/10^4, 'b =')"
+ ]
+ }
+],
+"metadata": {
+ "kernelspec": {
+ "display_name": "Scilab",
+ "language": "scilab",
+ "name": "scilab"
+ },
+ "language_info": {
+ "file_extension": ".sce",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "scilab",
+ "version": "0.7.1"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/6-Interpolation.ipynb b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/6-Interpolation.ipynb
new file mode 100644
index 0000000..9df0c59
--- /dev/null
+++ b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/6-Interpolation.ipynb
@@ -0,0 +1,1135 @@
+{
+"cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Chapter 6: Interpolation"
+ ]
+ },
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.10: Newton_Forward_Difference_Interpolation_Formula.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.10\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 1:5;\n",
+"Y = poly(0, 'Y');\n",
+"y = [2 5 7 Y 32];\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,5);\n",
+"del(:,1) = y';\n",
+"for j = 2:5\n",
+" for i = 1:n-j+1\n",
+" del(i,j) = del(i+1,j-1) - del(i,j-1);\n",
+" end\n",
+"end\n",
+"del(:,1) = [];\n",
+"\n",
+"// del4 = 0\n",
+"\n",
+"y0 = del(:,4);\n",
+"y0(isnan(y0)) = [];\n",
+"Y = roots(y0)\n",
+"disp(Y,'Missing value f(x3) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.11: Newton_Forward_Difference_Interpolation_Formula.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.11\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 0:5;\n",
+"y = [-3 3 11 27 57 107];\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,4);\n",
+"del(:,1) = y';\n",
+"for j = 2:4\n",
+" for i = 1:n-j+1\n",
+" del(i,j) = del(i+1,j-1) - del(i,j-1);\n",
+" end\n",
+"end\n",
+"del(:,1) = [];\n",
+"\n",
+"X = poly(0, 'x');\n",
+"h = x(2) - x(1);\n",
+"p = (X-x(1)) / h;\n",
+"x0 = x(1);\n",
+"y0 = y(1);\n",
+"dely0 = del(1,:);\n",
+"\n",
+"Y = y0;\n",
+"\n",
+"for i = 1:length(dely0)\n",
+" t = 1;\n",
+" for j = 1:i\n",
+" t = t * (p-j+1);\n",
+" end\n",
+" Y = Y + t*dely0(i)/factorial(i);\n",
+"end\n",
+"disp(Y,'Required cubic polynomial:')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.12: Newton_Backward_Difference_Interpolation_Formula.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.12\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 1:8;\n",
+"y = x^3;\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,4);\n",
+"del(:,1) = y';\n",
+"for j = 2:4\n",
+" for i = 1:n-j+1\n",
+" del(i+j-1,j) = del(i+j-1,j-1) - del(i+j-2,j-1);\n",
+" end\n",
+"end\n",
+"\n",
+"X = 7.5;\n",
+"h = x(2) - x(1);\n",
+"p = (X-x(n)) / h;\n",
+"xn = x(n);\n",
+"yn = y(n);\n",
+"delyn = del(n,:);\n",
+"\n",
+"Y = 0;\n",
+"\n",
+"for i = 0:length(delyn)-1\n",
+" t = 1;\n",
+" for j = 1:i\n",
+" t = t * (p+j-1);\n",
+" end\n",
+" Y = Y + t*delyn(i+1)/factorial(i);\n",
+"end\n",
+"disp(Y,'y(7.5) = ')\n",
+""
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.13: Newton_Backward_Difference_Interpolation_Formula.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.13\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 1974:2:1982;\n",
+"y = [40 43 48 52 57];\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,5);\n",
+"del(:,1) = y';\n",
+"for j = 2:5\n",
+" for i = 1:n-j+1\n",
+" del(i+j-1,j) = del(i+j-1,j-1) - del(i+j-2,j-1);\n",
+" end\n",
+"end\n",
+"\n",
+"X = 1979;\n",
+"h = x(2) - x(1);\n",
+"p = (X-x(n)) / h;\n",
+"xn = x(n);\n",
+"yn = y(n);\n",
+"delyn = del(n,:);\n",
+"\n",
+"Y = 0;\n",
+"\n",
+"for i = 0:length(delyn)-1\n",
+" t = 1;\n",
+" for j = 1:i\n",
+" t = t * (p+j-1);\n",
+" end\n",
+" Y = Y + t*delyn(i+1)/factorial(i);\n",
+"end\n",
+"Y = round(Y*10^4)/10^4;\n",
+"disp(Y,'Estimated sales for the year 1979: ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.14: Lagrange_Interpolation_Formula.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.14\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = [1 3 4 6];\n",
+"y = [-3 0 30 132];\n",
+"\n",
+"n = length(x);\n",
+"Y = 0;\n",
+"X = poly(0, 'X');\n",
+"//X = 5;\n",
+"for i = 1:n\n",
+" t = x;\n",
+" t(i) = [];\n",
+" p = 1;\n",
+" for j = 1:length(t)\n",
+" p = p * (X-t(j))/(x(i)-t(j));\n",
+" end\n",
+" Y = Y + p*y(i);\n",
+"end\n",
+"Y5 = horner(Y,5);\n",
+"disp(Y5,'y(5) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.15: Lagrange_Interpolation_Formula.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.15\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = [1 2 5];\n",
+"y = [1 4 10];\n",
+"\n",
+"n = length(x);\n",
+"Y = 0;\n",
+"X = poly(0, 'X');\n",
+"//X = 5;\n",
+"for i = 1:n\n",
+" t = x;\n",
+" t(i) = [];\n",
+" p = 1;\n",
+" for j = 1:length(t)\n",
+" p = p * (X-t(j))/(x(i)-t(j));\n",
+" end\n",
+" Y = Y + p*y(i);\n",
+"end\n",
+"Y5 = horner(Y,3);\n",
+"disp(Y5,'f(3) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.16: Lagrange_and_Newton_Divided_Difference_Interpolation_Formulae.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.16\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = [0 1 2 4];\n",
+"y = [1 1 2 5];\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,4);\n",
+"del(:,1) = y';\n",
+"for j = 2:4\n",
+" for i = 1:n-j+1\n",
+" del(i,j) = (del(i+1,j-1) - del(i,j-1)) / (x(i+j-1) - x(i));\n",
+" end\n",
+"end\n",
+"del(:,1) = [];\n",
+"\n",
+"Y = 0;\n",
+"X = poly(0, 'X');\n",
+"for i = 1:n\n",
+" t = x;\n",
+" t(i) = [];\n",
+" p = 1;\n",
+" for j = 1:length(t)\n",
+" p = p * (X-t(j))/(x(i)-t(j));\n",
+" end\n",
+" Y = Y + p*y(i);\n",
+"end\n",
+"disp(round(Y*10^4)/10^4,'Interpolating polynomial:')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.17: Newton_Divided_Difference_Interpolation_Formulae.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.17\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = [0 1 4];\n",
+"y = [2 1 4];\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,3);\n",
+"del(:,1) = y';\n",
+"for j = 2:3\n",
+" for i = 1:n-j+1\n",
+" del(i,j) = (del(i+1,j-1) - del(i,j-1)) / (x(i+j-1) - x(i));\n",
+" end\n",
+"end\n",
+"del(:,1) = [];\n",
+"\n",
+"Y = 0;\n",
+"X = 2;\n",
+"for i = 1:n\n",
+" t = x;\n",
+" t(i) = [];\n",
+" p = 1;\n",
+" for j = 1:length(t)\n",
+" p = p * (X-t(j))/(x(i)-t(j));\n",
+" end\n",
+" Y = Y + p*y(i);\n",
+"end\n",
+"disp(Y,'y(2) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.18: Identity_Proof_for_Newton_and_Lagrange_Interpolation_Formulae.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"// Example 6.18\n",
+"// This is an analytical problem and need not be coded."
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.19: Interpolation_in_Two_Dimensions.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.19\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 0:4;\n",
+"n = length(x);\n",
+"f = 'X^2 + Y^2 - Y';\n",
+"tab = %nan*ones(n,5);\n",
+"\n",
+"for j = 0:4\n",
+" fj = strsubst(f,'Y','j');\n",
+" for i = 1:n\n",
+" tab(i,j+1) = eval(strsubst(fj,'X','x(i)'));\n",
+" end\n",
+"end\n",
+"//tab(:,1) = [];\n",
+"mprintf('%4s %6s %6s %6s %6s %6s','x','y=0','y=1','y=2','y=3','y=4')\n",
+"disp([(0:4)' tab])\n",
+"tab2 = tab(2:4,2:4)';\n",
+"n1 = length(tab2(:,1));\n",
+"y = 2:4;\n",
+"\n",
+"del1 = %nan*ones(n1,3);\n",
+"del1(:,1) = tab2(:,1);\n",
+"for j = 2:4\n",
+" for i = 1:n1-j+1\n",
+" del1(i,j) = del1(i+1,j-1) - del1(i,j-1);\n",
+" end\n",
+"end\n",
+"\n",
+"del2 = %nan*ones(n1,3);\n",
+"del2(:,1) = tab2(:,2);\n",
+"for j = 2:4\n",
+" for i = 1:n1-j+1\n",
+" del2(i,j) = del2(i+1,j-1) - del2(i,j-1);\n",
+" end\n",
+"end\n",
+"\n",
+"del3 = %nan*ones(n1,3);\n",
+"del3(:,1) = tab2(:,3);\n",
+"for j = 2:4\n",
+" for i = 1:n1-j+1\n",
+" del3(i,j) = del3(i+1,j-1) - del3(i,j-1);\n",
+" end\n",
+"end\n",
+"\n",
+"y0 = y(1);\n",
+"Y = 3.5;\n",
+"hy = y(2) - y(1);\n",
+"py = (Y-y0)/hy;\n",
+"\n",
+"f1y = 0;\n",
+"del1y0 = del1(1,:);\n",
+"for i = 0:length(del1y0)-1\n",
+" t = 1;\n",
+" for j = 1:i\n",
+" t = t * (py-j+1);\n",
+" end\n",
+" f1y = f1y + t*del1y0(i+1)/factorial(i);\n",
+"end\n",
+"\n",
+"f2y = 0;\n",
+"del2y0 = del2(1,:);\n",
+"for i = 0:length(del2y0)-1\n",
+" t = 1;\n",
+" for j = 1:i\n",
+" t = t * (py-j+1);\n",
+" end\n",
+" f2y = f2y + t*del2y0(i+1)/factorial(i);\n",
+"end\n",
+"\n",
+"f3y = 0;\n",
+"del3y0 = del3(1,:);\n",
+"for i = 0:length(del3y0)-1\n",
+" t = 1;\n",
+" for j = 1:i\n",
+" t = t * (py-j+1);\n",
+" end\n",
+" f3y = f3y + t*del3y0(i+1)/factorial(i);\n",
+"end\n",
+"\n",
+"del = %nan*ones(n1,3);\n",
+"del(:,1) = [f1y; f2y; f3y];\n",
+"for j = 2:4\n",
+" for i = 1:n1-j+1\n",
+" del(i,j) = del(i+1,j-1) - del(i,j-1);\n",
+" end\n",
+"end\n",
+"\n",
+"f = 0;\n",
+"X = 2.5;\n",
+"x0 = x(2);\n",
+"hx = x(2) - x(1);\n",
+"px = (X-x0)/hx;\n",
+"del0 = del(1,:)\n",
+"for i = 0:length(del0)-1\n",
+" t = 1;\n",
+" for j = 1:i\n",
+" t = t * (px-j+1);\n",
+" end\n",
+" f = f + t*del0(i+1)/factorial(i);\n",
+"end\n",
+"disp(f,'f(2.5,3.5) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.1: Forward_Difference_Table.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.1\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 0.1:0.2:1.3;\n",
+"y = [0.003 0.067 0.148 0.248 0.37 0.518 0.697];\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,6);\n",
+"del(:,1) = y';\n",
+"for j = 2:6\n",
+" for i = 1:n-j+1\n",
+" del(i,j) = del(i+1,j-1) - del(i,j-1);\n",
+" end\n",
+"end\n",
+"del = [x' del];\n",
+"del = round(del*10^3)/10^3;\n",
+"mprintf('%5s %7s %8s %9s %8s %8s %8s','x','y','dy','d2y','d3y','d4y','d5y')\n",
+"disp(del)"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.20: Cubic_Spline_Curve.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.20\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [p] = cubicsplin(x,y)\n",
+"// Fits point data to cubic spline fit\n",
+"\n",
+"n = length(x);\n",
+"a = y(1:n-1); // Spline Initials\n",
+"\n",
+"M1 = zeros(3*(n-1));\n",
+"M2 = zeros(3*(n-1),1);\n",
+"// Point Substitutions\n",
+"for i = 1:n-1\n",
+" M1(i,i) = x(i+1) - x(i);\n",
+" M1(i,i+n-1) = (x(i+1) - x(i))^2;\n",
+" M1(i,i+2*(n-1)) = (x(i+1) - x(i))^3;\n",
+" M2(i) = y(i+1) - y(i);\n",
+"end\n",
+"\n",
+"// Knot equations\n",
+"for i = 1:n-2\n",
+" // Derivative (S') continuity\n",
+" M1(i+n-1,i) = 1;\n",
+" M1(i+n-1,i+1) = -1;\n",
+" M1(i+n-1,i+n-1) = 2*(x(i+1)-x(i));\n",
+" M1(i+n-1,i+2*(n-1)) = 3*(x(i+1)-x(i))^2;\n",
+" // S'' continuity\n",
+" M1(i+2*n-3,i+n-1) = 2;\n",
+" M1(i+2*n-3,i+n) = -2;\n",
+" M1(i+2*n-3,i+2*(n-1)) = 6*(x(i+1)-x(i));\n",
+"end\n",
+"// Given BC\n",
+"M1(3*n-4,n) = 1;\n",
+"M1(3*n-3,2*n-2) = 1;\n",
+"M1(3*n-3,3*n-3) = 3*(3-2);\n",
+"\n",
+"var = M1\M2;\n",
+"var = round(var);\n",
+"b = var(1:n-1);\n",
+"c = var(n:2*(n-1));\n",
+"d = var(2*(n-1)+1:3*(n-1));\n",
+"p = [d c b a(:)];\n",
+"endfunction\n",
+"\n",
+"x = 0:3;\n",
+"y = [1 4 0 -2];\n",
+"p = cubicsplin(x,y);\n",
+"for i = 1:length(p(:,1))\n",
+" disp(strcat(['S',string(i-1),'(x) =']))\n",
+" disp(poly(p(i,:),'X',['coeff']))\n",
+"end"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.21: Cubic_Spline_Curve.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.21\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [p] = cubicsplin(x,y)\n",
+"// Fits point data to cubic spline fit\n",
+"\n",
+"n = length(x);\n",
+"a = y(1:n-1); // Spline Initials\n",
+"\n",
+"M1 = zeros(3*(n-1));\n",
+"M2 = zeros(3*(n-1),1);\n",
+"// Point Substitutions\n",
+"for i = 1:n-1\n",
+" M1(i,i) = x(i+1) - x(i);\n",
+" M1(i,i+n-1) = (x(i+1) - x(i))^2;\n",
+" M1(i,i+2*(n-1)) = (x(i+1) - x(i))^3;\n",
+" M2(i) = y(i+1) - y(i);\n",
+"end\n",
+"\n",
+"// Knot equations\n",
+"for i = 1:n-2\n",
+" // Derivative (S') continuity\n",
+" M1(i+n-1,i) = 1;\n",
+" M1(i+n-1,i+1) = -1;\n",
+" M1(i+n-1,i+n-1) = 2*(x(i+1)-x(i));\n",
+" M1(i+n-1,i+2*(n-1)) = 3*(x(i+1)-x(i))^2;\n",
+" // S'' continuity\n",
+" M1(i+2*n-3,i+n-1) = 2;\n",
+" M1(i+2*n-3,i+n) = -2;\n",
+" M1(i+2*n-3,i+2*(n-1)) = 6*(x(i+1)-x(i));\n",
+"end\n",
+"// Given BC\n",
+"M1(3*n-4,1) = 1;\n",
+"M1(3*n-3,n-1) = 1;\n",
+"M1(3*n-3,2*n-2) = 2*(3-2);\n",
+"M1(3*n-3,3*n-3) = 3*(3-2)^2;\n",
+"M2(3*n-4) = 2;\n",
+"M2(3*n-3) = 2;\n",
+"\n",
+"var = M1\M2;\n",
+"var = round(var);\n",
+"b = var(1:n-1);\n",
+"c = var(n:2*(n-1));\n",
+"d = var(2*(n-1)+1:3*(n-1));\n",
+"\n",
+"p = [a(:) b c d];\n",
+"endfunction\n",
+"\n",
+"x = 0:3;\n",
+"y = [1 4 0 -2];\n",
+"p = cubicsplin(x,y);\n",
+"for i=1:length(p(:,1))\n",
+" disp(strcat(['S',string(i-1),'(x) = ']))\n",
+" disp(poly(p(i,:),'x',['coeff']))\n",
+"end"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.22: Minima_of_a_Tabulated_Function.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.22\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 3:8;\n",
+"y = [0.205 0.24 0.259 0.262 0.25 0.224];\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,5);\n",
+"del(:,1) = y';\n",
+"for j = 2:5\n",
+" for i = 1:n-j+1\n",
+" del(i,j) = del(i+1,j-1) - del(i,j-1);\n",
+" end\n",
+"end\n",
+"\n",
+"X = poly(0, 'X');\n",
+"x0 = x(1);\n",
+"y0 = y(1);\n",
+"h = x(2) - x(1);\n",
+"p = (X-x0)/h;\n",
+"del0 = del(1,:);\n",
+"del0 = round(del0*10^4)/10^4;\n",
+"del0 = del0(1:find(del0==0)-1);\n",
+"\n",
+"Y = 0;\n",
+"for i = 0:length(del0)-1\n",
+" t = 1;\n",
+" for j = 1:i\n",
+" t = t * (p-j+1);\n",
+" end\n",
+" Y = Y + t*del0(i+1)/factorial(i);\n",
+"end\n",
+"disp(Y,'y = ')\n",
+"\n",
+"dydx = derivat(Y);\n",
+"minx = roots(dydx);\n",
+"miny = round(horner(Y,minx)*10^5)/10^5;\n",
+"disp(minx,'min_x = ')\n",
+"disp(miny,'min_y = ')\n",
+"//min_y value is incorrectly displayed in textbook as 0.25425 instead of 0.26278"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.23: Maxima_of_a_Tabulated_Function.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.23\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = [-1 1 2 3];\n",
+"y = [-21 15 12 3];\n",
+"\n",
+"n = length(x);\n",
+"X = poly(0, 'X');\n",
+"Y = 0;\n",
+"for i = 1:n\n",
+" t = x;\n",
+" t(i) = [];\n",
+" p = 1;\n",
+" for j = 1:length(t)\n",
+" p = p * (X-t(j))/(x(i)-t(j));\n",
+" end\n",
+" Y = Y + p*y(i);\n",
+"end\n",
+"\n",
+"dydx = derivat(Y);\n",
+"extx = real(roots(dydx));\n",
+"extx = round(extx*10^4)/10^4;\n",
+"d2ydx = derivat(dydx);\n",
+"\n",
+"if horner(d2ydx,extx(1)) < 0 then\n",
+" maxx = extx(1);\n",
+" maxy = horner(Y,maxx);\n",
+"else\n",
+" maxx = extx(2);\n",
+" maxy = horner(Y,maxx);\n",
+"end\n",
+"maxy = round(maxy*10^4)/10^4;\n",
+"disp(maxx,'max_x = ')\n",
+"disp(maxy,'max_y = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.24: Determination_of_Function_Value.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.24\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 1:3:10;\n",
+"F = [500426 329240 175212 40365];\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,4);\n",
+"del(:,1) = F';\n",
+"for j = 2:4\n",
+" for i = 1:n-j+1\n",
+" del(i,j) = del(i+1,j-1) - del(i,j-1);\n",
+" end\n",
+"end\n",
+"\n",
+"del0 = del(1,:);\n",
+"X = 2;\n",
+"x0 = x(1);\n",
+"h = x(2) - x(1);\n",
+"p = (X-x0) / h;\n",
+"F2 = 0;\n",
+"for i = 0:length(del0)-1\n",
+" t = 1;\n",
+" for j = 1:i\n",
+" t = t * (p-j+1);\n",
+" end\n",
+" F2 = F2 + t*del0(i+1)/factorial(i);\n",
+"end\n",
+"\n",
+"f2 = F(1) - F2;\n",
+"disp(f2,'f(2) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.2: Expression_for_Finite_Difference_Elements.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"// Example 6.2\n",
+"// This is an analytical problem and need not be coded."
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.3: Expression_for_Finite_Difference_Elements.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"// Example 6.3\n",
+"// This is an analytical problem and need not be coded."
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.4: Expression_for_Finite_Difference_Elements.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"// Example 6.4\n",
+"// This is an analytical problem and need not be coded."
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.5: Proof_of_Relatio.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"// Example 6.5\n",
+"// This is an analytical problem and need not be coded."
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.6: Proofs_of_given_Relations.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"// Example 6.6\n",
+"// This is an analytical problem and need not be coded."
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.7: Proof_for_Commutation_of_given_Operations.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"// Example 6.7\n",
+"// This is an analytical problem and need not be coded."
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.8: Newton_Forward_Difference_Interpolation_Formula.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.8\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 10:10:50;\n",
+"y = [46 66 81 93 101];\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,5);\n",
+"del(:,1) = y';\n",
+"for j = 2:5\n",
+" for i = 1:n-j+1\n",
+" del(i,j) = del(i+1,j-1) - del(i,j-1);\n",
+" end\n",
+"end\n",
+"del(:,1) = [];\n",
+"\n",
+"X = 15; //input\n",
+"for i = 1:n\n",
+" if X>x(i) then\n",
+" h = x(i+1) - x(i);\n",
+" p = (X-x(i)) / h;\n",
+" x0 = x(i);\n",
+" y0 = y(i);\n",
+" dely0 = del(i,:);\n",
+" dely0(isnan(y0)) = [];\n",
+" end\n",
+"end\n",
+"\n",
+"Y = y0;\n",
+"\n",
+"for i = 1:length(dely0)\n",
+" t = 1;\n",
+" for j = 1:i\n",
+" t = t * (p-j+1);\n",
+" end\n",
+" Y = Y + t*dely0(i)/factorial(i);\n",
+"end\n",
+"Y = round(Y*10^4)/10^4;\n",
+"disp(Y,'f(15) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 6.9: Newton_Forward_Difference_Interpolation_Formula.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 6.9\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 0.1:0.1:0.5;\n",
+"y = [1.4 1.56 1.76 2 2.28];\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,5);\n",
+"del(:,1) = y';\n",
+"for j = 2:5\n",
+" for i = 1:n-j+1\n",
+" del(i,j) = del(i+1,j-1) - del(i,j-1);\n",
+" end\n",
+"end\n",
+"del(:,1) = [];\n",
+"\n",
+"X = poly(0, 'X');\n",
+"h = x(2) - x(1);\n",
+"p = (X-x(1)) / h;\n",
+"x0 = x(1);\n",
+"y0 = y(1);\n",
+"dely0 = del(1,:);\n",
+"\n",
+"Y = y0;\n",
+"\n",
+"for i = 1:length(dely0)\n",
+" t = 1;\n",
+" for j = 1:i\n",
+" t = t * (p-j+1);\n",
+" end\n",
+" Y = Y + t*dely0(i)/factorial(i);\n",
+"end\n",
+"Y = round(Y*10^2)/10^2;\n",
+"disp(Y,'Required Newton''s Interpolating Polynomial:')"
+ ]
+ }
+],
+"metadata": {
+ "kernelspec": {
+ "display_name": "Scilab",
+ "language": "scilab",
+ "name": "scilab"
+ },
+ "language_info": {
+ "file_extension": ".sce",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "scilab",
+ "version": "0.7.1"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/7-Numerical_Differentiation_and_Integration.ipynb b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/7-Numerical_Differentiation_and_Integration.ipynb
new file mode 100644
index 0000000..3ad5aff
--- /dev/null
+++ b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/7-Numerical_Differentiation_and_Integration.ipynb
@@ -0,0 +1,887 @@
+{
+"cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Chapter 7: Numerical Differentiation and Integration"
+ ]
+ },
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.10: Integral_using_Simpson_One_Third_Rule.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Etample 7.10\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"t = 0:10:80;\n",
+"a = [30 31.63 33.34 35.47 37.75 40.33 43.25 46.69 50.67];\n",
+"\n",
+"h = t(2) - t(1);\n",
+"n = length(t);\n",
+"\n",
+"Is13 = a(1);\n",
+"for i = 2:n-1\n",
+" rem2 = i-fix(i./2).*2;\n",
+" if rem2 == 0 then\n",
+" Is13 = Is13 + 4*a(i);\n",
+" else\n",
+" Is13 = Is13 + 2*a(i);\n",
+" end\n",
+"end\n",
+"Is13 = (Is13 + a(n))/10^3;\n",
+"Is13 = round(h/3*Is13*10^4)/10^4;\n",
+"disp(strcat(['v = ',string(Is13),' km/s']))"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.11: Romberg_Integration_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.11\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 1:0.1:1.8;\n",
+"x = round(x*10)/10;\n",
+"y = [1.543 1.669 1.811 1.971 2.151 2.352 2.577 2.828 3.107];\n",
+"n = length(x);\n",
+"x0 = x(1);\n",
+"xn = x(n);\n",
+"\n",
+"N = [1 2 4 8]\n",
+"for j = 1:length(N)\n",
+" h = (xn - x0)./N(j);\n",
+" I = y(1);\n",
+" for xx = x0+h:h:xn-h\n",
+" xx = round(xx*10)/10;\n",
+" I = I + 2*y(x==xx);\n",
+" end\n",
+" Itrap(j) = h/2*(I + y(n));\n",
+" IRomb(1) = Itrap(1);\n",
+" if j~=1 then\n",
+" IRomb(j) = (4^(j-1)*Itrap(j)-IRomb(j-1))/(4^(j-1)-1);\n",
+" end\n",
+"end\n",
+"IRomb = round(IRomb*10^5)/10^5;\n",
+"\n",
+"disp(Itrap(length(N)),'Integral using Trapezoidal rule:')\n",
+"disp(IRomb(length(N)),'Integral using Romberg''s formula:')\n",
+"//In third step of computation of integral using Romberg's formula, author mistakenly took the 1.7672 instead of 1.7684 which resulted in a difference"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.12: Double_Integral_using_Trapezoidal_Rule.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.12\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f] = fun1(x,y)\n",
+" f = 1 / (x+y);\n",
+"endfunction\n",
+"\n",
+"x = 1:0.25:2;\n",
+"y = x;\n",
+"\n",
+"m = length(x);\n",
+"n = length(y);\n",
+"\n",
+"del = %nan*ones(m,n);\n",
+"for j = 1:n\n",
+" for i = 1:m\n",
+" del(i,j) = fun1(x(i),y(j));\n",
+" end\n",
+"end\n",
+"\n",
+"hx = x(2) - x(1);\n",
+"for i = 1:m\n",
+" I = del(i,1);\n",
+" for j = 2:n-1\n",
+" I = I + 2*del(i,j);\n",
+" end\n",
+" Itrap1(i) = hx/2 * (I+del(i,n));\n",
+"end\n",
+"Itrap1 = round(Itrap1*10^4)/10^4;\n",
+"\n",
+"hy = y(2) - y(1);\n",
+"Itrap2 = Itrap1(1)\n",
+"for i = 2:n-1\n",
+" Itrap2 = Itrap2 + 2* Itrap1(i);\n",
+"end\n",
+"Itrap2 = round(hy/2*(Itrap2+Itrap1(m))*10^4)/10^4;\n",
+"disp(Itrap2,'I = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.13: Double_Integral_using_Trapezoidal_Rule.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.13\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f] = fun1(x,y)\n",
+" f = sqrt(sin(x+y));\n",
+"endfunction\n",
+"\n",
+"x = 0:%pi/8:%pi/2;\n",
+"y = x;\n",
+"\n",
+"m = length(x);\n",
+"n = length(y);\n",
+"\n",
+"del = %nan*ones(m,n);\n",
+"for j = 1:n\n",
+" for i = 1:m\n",
+" del(i,j) = fun1(x(i),y(j));\n",
+" end\n",
+"end\n",
+"\n",
+"hx = x(2) - x(1);\n",
+"for i = 1:m\n",
+" I = del(i,1);\n",
+" for j = 2:n-1\n",
+" I = I + 2*del(i,j);\n",
+" end\n",
+" Itrap1(i) = hx/2 * (I+del(i,n));\n",
+"end\n",
+"Itrap1 = round(Itrap1*10^4)/10^4;\n",
+"\n",
+"hy = y(2) - y(1);\n",
+"Itrap2 = Itrap1(1)\n",
+"for i = 2:n-1\n",
+" Itrap2 = Itrap2 + 2* Itrap1(i);\n",
+"end\n",
+"Itrap2 = round(hy/2*(Itrap2+Itrap1(m))*10^4)/10^4;\n",
+"disp(Itrap2,'I = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.14: One_Point_Gauss_Legendre_Quadrature_Formula.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.14\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"n = 1;\n",
+"if n==1 then\n",
+" M = [0 2];\n",
+"elseif n==2\n",
+" M = [sqrt(1/3) 1; -sqrt(1/3) 1];\n",
+"elseif n==3\n",
+" M = [0 8/9; -0.774597 5/9; 0.774597 5/9];\n",
+"elseif n==4\n",
+" M = [-0.339981 0.652145; -0.861136 0.347855; 0.339981 0.652145; 0.861136 0.347855];\n",
+"elseif n==5\n",
+" M = [-0 0.568889; -0.538469 0.467914; -0.906180 0.236927; 0 0.568889; 0.538469 0.467914; 0.906180 0.236927];\n",
+"elseif n==6\n",
+" M = [-0.238619 0.467914; -0.661209 0.360762; -0.932470 0.171325; 0.238619 0.467914; 0.661209 0.360762; 0.932470 0.171325];\n",
+"end\n",
+"\n",
+"X = M(:,1);\n",
+"W = M(:,2);\n",
+"\n",
+"disp(X,'E1 = ')\n",
+"disp(W,'W1 = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.15: Two_Point_Gauss_Legendre_Quadrature_Formula.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.15\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"n = 2;\n",
+"if n==1 then\n",
+" M = [0 2];\n",
+"elseif n==2\n",
+" M = [sqrt(1/3) 1; -sqrt(1/3) 1];\n",
+"elseif n==3\n",
+" M = [0 8/9; -0.774597 5/9; 0.774597 5/9];\n",
+"elseif n==4\n",
+" M = [-0.339981 0.652145; -0.861136 0.347855; 0.339981 0.652145; 0.861136 0.347855];\n",
+"elseif n==5\n",
+" M = [-0 0.568889; -0.538469 0.467914; -0.906180 0.236927; 0 0.568889; 0.538469 0.467914; 0.906180 0.236927];\n",
+"elseif n==6\n",
+" M = [-0.238619 0.467914; -0.661209 0.360762; -0.932470 0.171325; 0.238619 0.467914; 0.661209 0.360762; 0.932470 0.171325];\n",
+"end\n",
+"\n",
+"X = M(:,1);\n",
+"W = M(:,2);\n",
+"\n",
+"disp(W(1),'W1 = ')\n",
+"disp(W(2),'W2 = ')\n",
+"disp(X(1),'E1 = ')\n",
+"disp(X(2),'E2 = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.16: Four_Point_Gauss_Legendre_Quadrature_Formula.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.16\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f] = fun1(x)\n",
+" f = 3*x^2 + x^3;\n",
+"endfunction\n",
+"\n",
+"n = 4;\n",
+"if n==1 then\n",
+" M = [0 2];\n",
+"elseif n==2\n",
+" M = [sqrt(1/3) 1];\n",
+"elseif n==3\n",
+" M = [0 8/9; -0.774597 5/9; 0.774597 5/9];\n",
+"elseif n==4\n",
+" M = [-0.339981 0.652145; -0.861136 0.347855; 0.339981 0.652145; 0.861136 0.347855];\n",
+"elseif n==5\n",
+" M = [-0 0.568889; -0.538469 0.467914; -0.906180 0.236927; 0 0.568889; 0.538469 0.467914; 0.906180 0.236927];\n",
+"elseif n==6\n",
+" M = [-0.238619 0.467914; -0.661209 0.360762; -0.932470 0.171325; 0.238619 0.467914; 0.661209 0.360762; 0.932470 0.171325];\n",
+"end\n",
+"\n",
+"X = M(:,1);\n",
+"W = M(:,2);\n",
+"I = 0;\n",
+"for i = 1:length(X)\n",
+" I = I + W(i)*fun1(X(i));\n",
+"end\n",
+"disp(I,'I = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.1: Determination_of_Differential_Function_Value.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.1\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 0:0.2:1;\n",
+"y = [1 1.16 3.56 13.96 41.96 101];\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,6);\n",
+"del(:,1) = y';\n",
+"for j = 2:6\n",
+" for i = 1:n-j+1\n",
+" del(i,j) = del(i+1,j-1) - del(i,j-1);\n",
+" end\n",
+"end\n",
+"del = round(del*10^2)/10^2;\n",
+"mprintf('%5s %6s %9s %8s %8s %8s %7s','x','y','dy','d2y','d3y','d4y','d5y')\n",
+"disp([x' del])\n",
+"\n",
+"h = x(2) - x(1);\n",
+"del0 = del(1,:);\n",
+"del1 = del(2,:);\n",
+"\n",
+"df1 = (del1(2) - del1(3)/2 + del1(4)/3 - del1(5)/4) / h;\n",
+"d2f0 = (del0(2) - del0(3) + del0(4)*11/12 - del0(5)*5/6) / h^2;\n",
+"disp(round(d2f0*10^1)/10^1,'f''''(0) = ')\n",
+"disp(round(df1*10)/10,'f''(0.2) = ')\n",
+""
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.2: Determination_of_Differential_Function_Value.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.2\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 1.4:0.2:2.2;\n",
+"y = [4.0552 4.953 6.0496 7.3891 9.025];\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,5);\n",
+"del(:,1) = y';\n",
+"for j = 2:5\n",
+" for i = 1:n-j+1\n",
+" del(i+j-1,j) = del(i+j-1,j-1) - del(i+j-2,j-1);\n",
+" end\n",
+"end\n",
+"mprintf('%5s %6s %10s %10s %9s %9s','x','y','dy','d2y','d3y','d4y')\n",
+"disp([x' del])\n",
+"\n",
+"h = x(2) - x(1);\n",
+"deln = del(5,:);\n",
+"\n",
+"dfn = (deln(2) + deln(3)/2 + deln(4)/3 + deln(5)/4) / h;\n",
+"d2fn = (deln(3) + deln(4) + deln(5)*11/12) / h^2;\n",
+"dfn = round(dfn*10^4)/10^4;\n",
+"d2fn = round(d2fn*10^4)/10^4;\n",
+"disp(dfn,'y''(2.2) = ')\n",
+"disp(d2fn,'y''''(2.2) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.3: Determination_of_Differential_Function_Value.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.3\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = 0:4;\n",
+"y = [6.9897 7.4036 7.7815 8.1281 8.451];\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,5);\n",
+"del(:,1) = y';\n",
+"for j = 2:6\n",
+" for i = 1:n-j+1\n",
+" del(i,j) = del(i+1,j-1) - del(i,j-1);\n",
+" end\n",
+"end\n",
+"del(:,1) = [];\n",
+"n0 = length(del(1,:));\n",
+"\n",
+"X = 2;\n",
+"i = find(x==X);\n",
+"dowy = 0;\n",
+"\n",
+"for j = 1:n0\n",
+" if j==2*int(j/2) then\n",
+" add = del(i,j);\n",
+" else\n",
+" add = (del(i-1,j) + del(i,j))/2;\n",
+" i = i-1;\n",
+" if i==0 then\n",
+" break\n",
+" end\n",
+" end\n",
+"\n",
+" if add == %nan then\n",
+" break\n",
+" else\n",
+" dowy(j) = add;\n",
+" end\n",
+"end\n",
+"mprintf('%5s %6s %10s %9s %9s %9s','x','y','dy','d2y','d3y','d4y')\n",
+"disp([x' y' del])\n",
+"\n",
+"mu = 1;\n",
+"h = x(2) - x(1);\n",
+"dy2 = mu/h*(dowy(1) - 1/6*dowy(3));\n",
+"d2y2 = mu/h^2*(dowy(2)-1/12*dowy(4));\n",
+"dy2 = round(dy2*10^4)/10^4;\n",
+"d2y2 = round(d2y2*10^4)/10^4;\n",
+"\n",
+"disp(dy2,'y''(2) = ')\n",
+"disp(d2y2,'y''''(2) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.4: Determination_of_Differential_Function_Value.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.4\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"x = [0.15 0.21 0.23 0.27 0.32 0.35];\n",
+"y = [0.1761 0.3222 0.3617 0.4314 0.5051 0.5441];\n",
+"\n",
+"n = length(x);\n",
+"del = %nan*ones(n,6);\n",
+"del(:,1) = y';\n",
+"for j = 2:6\n",
+" for i = 1:n-j+1\n",
+" del(i,j) = (del(i+1,j-1) - del(i,j-1)) / (x(i+j-1)-x(i));\n",
+" end\n",
+"end\n",
+"del(:,1) = [];\n",
+"del = round(del*10^3)/10^3;\n",
+"mprintf('%5s %6s %10s %10s %8s %9s %9s','x','y','dy','d2y','d3y','d4y','d5y')\n",
+"disp([x' y' del])\n",
+"X = poly(0, 'X');\n",
+"del0 = del(1,:);\n",
+"y0 = y(1);\n",
+"Y = y0;\n",
+"for i = 1:length(del0)\n",
+" p = 1;\n",
+" for j = 1:i\n",
+" p = p*(X-x(j));\n",
+" end\n",
+" Y = Y + p*del0(i);\n",
+"end\n",
+"\n",
+"dydx = derivat(Y);\n",
+"d2ydx = derivat(dydx);\n",
+"\n",
+"XX = 0.25;\n",
+"dy = horner(dydx,XX);\n",
+"d2y = horner(d2ydx,XX);\n",
+"\n",
+"disp(round(dy*10^4)/10^4,'y''(0.25) = ')\n",
+"disp(d2ydx,'y''''(x) = ')\n",
+"disp(d2y,'y''''(0.25) = ')\n",
+"//The constant term in y''(x) is incorrectly computed to -91.7 instead of -97.42 in the text."
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.5: Richardson_Extrapolation_Limit.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.5\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f] = y(x)\n",
+" f = -1/x;\n",
+"endfunction\n",
+"\n",
+"H = [0.0128 0.0064 0.0032];\n",
+"n = length(H);\n",
+"x = 0.05;\n",
+"h = H(1);\n",
+"Fh = (y(x+h) - y(x-h)) / (2*h);\n",
+"Fh2 = (y(x+h/2) - y(x-h/2)) / (h);\n",
+"Fh4 = (y(x+h/4) - y(x-h/4)) / (h/2);\n",
+"\n",
+"F1h2 = (4*Fh2 - Fh) / (4-1);\n",
+"F1h4 = (4*Fh4 - Fh2) / (4-1);\n",
+"F2h4 = (4^2*F1h4 - F1h2) / (4^2-1);\n",
+"del = %nan*ones(n,3);\n",
+"del(:,1) = [Fh Fh2 Fh4]';\n",
+"del(1:2,2) = [F1h2 F1h4]';\n",
+"del(1,3) = F2h4;\n",
+"\n",
+"disp(del(1,n),'y''(0.05) = ')\n",
+"Exact = 1/x^2;\n",
+"disp(Exact,'Exact Value:')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.6: Integral_using_Trapezoidal_and_Simpson_One_Third_Rule.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.6\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [I] = trap (fun,a,b,n)\n",
+"// Integrate the function over the interval using Trapezoidal Formula\n",
+"// trap (fun,a,b,n)\n",
+"// fun - function to be integrated\n",
+"// a - lower limit of integration\n",
+"// b - upper limit of integration\n",
+"// n - No. of times trapezoidal rule needs to be performed\n",
+"\n",
+"N = n + 1; // N - total no. of points\n",
+"h = (b-a) / (N-1);\n",
+"x = linspace(a,b,N);\n",
+"y = fun(x);\n",
+"\n",
+"sum1 = y(1) + 2 * sum(y(2:N-1)) + y(N);\n",
+"I = h * sum1 / 2; // Trapezoidal Integral Value\n",
+"endfunction\n",
+"\n",
+"function [I] = simp13 (fun,a,b,n)\n",
+"// Integrate the function over the interval using Simpson's 1/3rd rule\n",
+"// simp13 (fun,a,b,n)\n",
+"// fun - function to be integrated\n",
+"// a - lower limit of integration\n",
+"// b - upper limit of integration\n",
+"// n - No. of times simpson's 1/3rd rule needs to be performed\n",
+"\n",
+"N = 2 * n + 1; // N - total no. of points\n",
+"h = (b-a) / (N-1);\n",
+"x = linspace(a,b,N);\n",
+"y = fun(x);\n",
+"\n",
+"sum1 = y(1) + 4 * sum(y(2:2:N-1)) + 2 * sum(y(3:2:N-2)) + y(N);\n",
+"I = h* sum1 / 3; // Simpson's 1/3rd Integral Value\n",
+"endfunction\n",
+"\n",
+"n = 6;\n",
+"ntrap = n;\n",
+"ns13 = n/2;\n",
+"I = [trap(sin,0,%pi,ntrap); simp13(sin,0,%pi,ns13)];\n",
+"I = round(I*10^4)/10^4;\n",
+"true = integrate('sin(x)','x',0,%pi);\n",
+"err = abs(true - I) / true*100;\n",
+"err = round(err*100)/100;\n",
+"\n",
+"disp(I(1),'y_trap = ')\n",
+"disp(I(2),'y_simp13 = ')\n",
+"disp(err(1),'error_trap = ')\n",
+"disp(err(2),'error_simp13 = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.7: Integral_using_Simpson_One_Third_Rule.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.7\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [I] = simp13 (fun,a,b,n)\n",
+"// Integrate the function over the interval using Simpson's 1/3rd rule\n",
+"// simp13 (fun,a,b,n)\n",
+"// fun - function to be integrated\n",
+"// a - lower limit of integration\n",
+"// b - upper limit of integration\n",
+"// n - No. of times simpson's 1/3rd rule needs to be performed\n",
+"\n",
+"N = 2 * n + 1; // N - total no. of points\n",
+"h = (b-a) / (N-1);\n",
+"x = linspace(a,b,N);\n",
+"y = fun(x);\n",
+"\n",
+"sum1 = y(1) + 4 * sum(y(2:2:N-1)) + 2 * sum(y(3:2:N-2)) + y(N);\n",
+"I = h* sum1 / 3; // Simpson's 1/3rd Integral Value\n",
+"endfunction\n",
+"\n",
+"n = 8;\n",
+"ns13 = n/2;\n",
+"I = simp13(log,1,5,ns13);\n",
+"I = round(I*10^4)/10^4;\n",
+"deff('[y] = true(x)',['y = x * log(x) - x']);\n",
+"trueVal = true(5) - true(1);\n",
+"err = abs(trueVal - I) / trueVal*100;\n",
+"err = round(err*100)/100;\n",
+"\n",
+"disp(I,'y_simp13 = ')\n",
+"disp(trueVal,'Actual Integral = ')\n",
+"disp(err,'error_simp13 = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.8: Integral_using_Trapezoidal_and_Simpson_One_Third_Rule.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.8\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [I] = trap (fun,a,b,n)\n",
+"// Integrate the function over the interval using Trapezoidal Formula\n",
+"// trap (fun,a,b,n)\n",
+"// fun - function to be integrated\n",
+"// a - lower limit of integration\n",
+"// b - upper limit of integration\n",
+"// n - No. of times trapezoidal rule needs to be performed\n",
+"\n",
+"N = n + 1; // N - total no. of points\n",
+"h = (b-a) / (N-1);\n",
+"x = linspace(a,b,N);\n",
+"y = fun(x);\n",
+"\n",
+"sum1 = y(1) + 2 * sum(y(2:N-1)) + y(N);\n",
+"I = h * sum1 / 2; // Trapezoidal Integral Value\n",
+"endfunction\n",
+"\n",
+"function [I] = simp13 (fun,a,b,n)\n",
+"// Integrate the function over the interval using Simpson's 1/3rd rule\n",
+"// simp13 (fun,a,b,n)\n",
+"// fun - function to be integrated\n",
+"// a - lower limit of integration\n",
+"// b - upper limit of integration\n",
+"// n - No. of times simpson's 1/3rd rule needs to be performed\n",
+"\n",
+"N = 2 * n + 1; // N - total no. of points\n",
+"h = (b-a) / (N-1);\n",
+"x = linspace(a,b,N);\n",
+"y = fun(x);\n",
+"\n",
+"sum1 = y(1) + 4 * sum(y(2:2:N-1)) + 2 * sum(y(3:2:N-2)) + y(N);\n",
+"I = h* sum1 / 3; // Simpson's 1/3rd Integral Value\n",
+"endfunction\n",
+"\n",
+"function [f] = fun1(x)\n",
+" f = 1 ./(1+x^2);\n",
+"endfunction\n",
+"\n",
+"\n",
+"n = 4;\n",
+"ntrap = n;\n",
+"ns13 = n/2;\n",
+"I = [trap(fun1,0,1,ntrap); simp13(fun1,0,1,ns13)];\n",
+"I = round(I*10^4)/10^4;\n",
+"true = intg(0,1,fun1);\n",
+"\n",
+"disp(I(1),'y_trap = ')\n",
+"disp(I(2),'y_simp13 = ')\n",
+"disp(I(2)*4,'Approx pi = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 7.9: Integral_using_Simpson_One_Third_Rule.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 7.9\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [I] = simp13 (fun,a,b,n)\n",
+"// Integrate the function over the interval using Simpson's 1/3rd rule\n",
+"// simp13 (fun,a,b,n)\n",
+"// fun - function to be integrated\n",
+"// a - lower limit of integration\n",
+"// b - upper limit of integration\n",
+"// n - No. of times simpson's 1/3rd rule needs to be performed\n",
+"\n",
+"N = 2 * n + 1; // N - total no. of points\n",
+"h = (b-a) / (N-1);\n",
+"x = linspace(a,b,N);\n",
+"y = fun(x);\n",
+"\n",
+"sum1 = y(1) + 4 * sum(y(2:2:N-1)) + 2 * sum(y(3:2:N-2)) + y(N);\n",
+"I = h* sum1 / 3; // Simpson's 1/3rd Integral Value\n",
+"endfunction\n",
+"\n",
+"function [f] = fun1(x)\n",
+" f = sqrt(2/%pi)*exp(-x^2/2);\n",
+"endfunction\n",
+"\n",
+"h = 0.125;\n",
+"n = (1-0)/h;\n",
+"ns13 = n/2;\n",
+"I = simp13(fun1,0,1,ns13);\n",
+"I = round(I*10^4)/10^4;\n",
+"\n",
+"disp(I,'Integral value, I = ')"
+ ]
+ }
+],
+"metadata": {
+ "kernelspec": {
+ "display_name": "Scilab",
+ "language": "scilab",
+ "name": "scilab"
+ },
+ "language_info": {
+ "file_extension": ".sce",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "scilab",
+ "version": "0.7.1"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/8-Ordinary_Differential_Equations.ipynb b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/8-Ordinary_Differential_Equations.ipynb
new file mode 100644
index 0000000..9abaf29
--- /dev/null
+++ b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/8-Ordinary_Differential_Equations.ipynb
@@ -0,0 +1,497 @@
+{
+"cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Chapter 8: Ordinary Differential Equations"
+ ]
+ },
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 8.1: Initial_Value_Problem_using_Taylor_Series_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 8.1\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f] = dydt(t,y)\n",
+" f = t+y;\n",
+"endfunction\n",
+"\n",
+"y0 = 0;\n",
+"t0 = 1;\n",
+"t = 1.2;\n",
+"h = 0.1;\n",
+"\n",
+"n = (t-t0)/h;\n",
+"tt = t0;\n",
+"y = y0;\n",
+"den = [1 2 6 24 120];\n",
+"for i = 1:n\n",
+" d2ydt = 1 + dydt(tt,y);\n",
+" d3ydt = d2ydt;\n",
+" d4ydt = d3ydt;\n",
+" d5ydt = d4ydt;\n",
+" dy = [dydt(tt,y) d2ydt d3ydt d4ydt d5ydt];\n",
+" tt = tt + h;\n",
+" for j = 1:length(dy)\n",
+" y = y + dy(j)*(tt-t0)^j/den(j);\n",
+" end\n",
+" t0 = tt;\n",
+"end\n",
+"disp(y,'y(1.2) = ')\n",
+"\n",
+"function [f] = closed(t)\n",
+" f = -t -1 + 2*exp(t-1);\n",
+"endfunction\n",
+"yclosed = closed(1.2);\n",
+"yclosed = round(yclosed*10^4)/10^4;\n",
+"disp(yclosed,'y_closed form = ')\n",
+"disp('Comparing the results obtained numerically and in closed form, we observe ')\n",
+"disp('that they agree up to four decimals')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 8.2: Initial_Value_Problem_using_Euler_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 8.2\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f] = dydt(t,y)\n",
+" f = (y-t) / (y+t);\n",
+"endfunction\n",
+"\n",
+"y0 = 1;\n",
+"t0 = 0;\n",
+"t = 0.1;\n",
+"n = 5;\n",
+"h = (t-t0)/n;\n",
+"\n",
+"tt = t0;\n",
+"y = y0;\n",
+"for i = 1:n\n",
+" y = y +h*dydt(tt,y);\n",
+" y = round(y*10^4)/10^4;\n",
+" tt = tt + h;\n",
+"end\n",
+"disp(y,'y(t = 0.1) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 8.3: Initial_Value_Problem_using_Modified_Euler_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 8.3\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f] = dydt(t,y)\n",
+" f = t + sqrt(y);\n",
+"endfunction\n",
+"\n",
+"y0 = 1;\n",
+"t0 = 0;\n",
+"h = 0.2;\n",
+"t = 0.6;\n",
+"n = (t-t0)/h;\n",
+"\n",
+"tt = t0;\n",
+"\n",
+"for i = 1:n\n",
+" y11 = y0 + h*dydt(tt,y0);\n",
+" t1 = tt + h;\n",
+" y1 = y0 + h/2*(dydt(tt,y0) + dydt(t1,y11));\n",
+" y1 = round(y1*10^4)/10^4;\n",
+"\n",
+" y(i) = y1;\n",
+" y0 = y1;\n",
+" tt = t1;\n",
+"end\n",
+"mprintf('%5s %8s','t','y')\n",
+"disp([(t0+h:h:t)' y])"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 8.4: Initial_Value_Problem_using_Second_Order_Runge_Kutta_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 8.4\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f] = fun1(x,y)\n",
+" f = (y+x) / (y-x);\n",
+"endfunction\n",
+"\n",
+"function [f] = rk2(x,y)\n",
+" k1 = h*fun1(x,y);\n",
+" k2 = h*fun1(x+3/2*h,y+3/2*k1);\n",
+" f = y + 1/3*(2*k1+k2);\n",
+"endfunction\n",
+"\n",
+"x0 = 0;\n",
+"y0 = 1;\n",
+"h = 0.2;\n",
+"x = 0.4;\n",
+"n = (x-x0)/h;\n",
+"\n",
+"for i = 1:n\n",
+" y = rk2(x0,y0);\n",
+" x0 = x0 + h;\n",
+" y0 = y;\n",
+" y = round(y*10^5)/10^5;\n",
+"end\n",
+"\n",
+"disp(y,'y(0.4) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 8.5: Initial_Value_Problem_using_Fourth_Order_Runge_Kutta_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Etample 8.5\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f] = fun1(t,y)\n",
+" f = t+y;\n",
+"endfunction\n",
+"\n",
+"function [f] = rk4(t,y)\n",
+" k1 = h*fun1(t,y);\n",
+" k2 = h*fun1(t+1/2*h,y+1/2*k1);\n",
+" k3 = h*fun1(t+1/2*h,y+1/2*k2);\n",
+" k4 = h*fun1(t+h,y+k1);\n",
+" f = y + 1/6*(k1+2*k2+2*k3+k4);\n",
+"endfunction\n",
+"\n",
+"t0 = 0;\n",
+"y0 = 1;\n",
+"h = 0.1;\n",
+"t = 0.4;\n",
+"n = (t-t0)/h;\n",
+"\n",
+"for i = 1:n\n",
+" y = rk4(t0,y0);\n",
+" t0 = t0 + h;\n",
+" y0 = y;\n",
+" y = round(y*10^5)/10^5;\n",
+"end\n",
+"\n",
+"disp(y,'y(0.4) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 8.6: Van_Der_Pol_Equation_using_Fourth_Order_Runge_Kutta_Equation.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 8.6\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f] = f1(x,y,p)\n",
+" f = p;\n",
+"endfunction\n",
+"\n",
+"function [f] = f2(x,y,p)\n",
+" f = 0.1*(1-y^2)*p - y;\n",
+"endfunction\n",
+"\n",
+"x0 = 0;\n",
+"y0 = 1;\n",
+"p0 = 0;\n",
+"h = 0.2;\n",
+"x = 0.2;\n",
+"n = (x-x0)/h;\n",
+"\n",
+"for i = 1:n\n",
+" k1 = h*f1(x0,y0,p0);\n",
+" l1 = h*f2(x0,y0,p0);\n",
+" k2 = h*f1(x0+h/2,y0+k1/2,p0+l1/2);\n",
+" l2 = h*f2(x0+h/2,y0+k1/2,p0+l1/2);\n",
+" k3 = h*f1(x0+h/2,y0+k2/2,p0+l2/2);\n",
+" l3 = h*f2(x0+h/2,y0+k2/2,p0+l2/2);\n",
+" k4 = h*f1(x0+h,y0+k3,p0+l3);\n",
+" l4 = h*f2(x0+h,y0+k3,p0+l3);\n",
+" y = y0 + 1/6*(k1+2*(k2+k3)+k4);\n",
+" p = p0 + 1/6*(l1+2*(l2+l3)+l4);\n",
+" y = round(y*10^4)/10^4;\n",
+" p = round(p*10^4)/10^4;\n",
+"end\n",
+"\n",
+"disp(y,'y(0.2) = ')\n",
+"disp(p,'y''(0.2) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 8.7: Milne_Predictor_Corrector_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 8.7\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f] = dy(t,y)\n",
+" f = 1/2*(t+y);\n",
+"endfunction\n",
+"\n",
+"tt = 0:0.5:1.5;\n",
+"yy = [2 2.636 3.595 4.968];\n",
+"\n",
+"t0 = tt(1);\n",
+"y0 = yy(1);\n",
+"t = 2;\n",
+"h = tt(2) - tt(1);\n",
+"n = (t-t0)/h;\n",
+"for i = 1:n\n",
+" dydt(1) = dy(t0,yy(1));\n",
+" dydt(2) = dy(t0+h,yy(2));\n",
+" dydt(3) = dy(t0+2*h,yy(3));\n",
+" dydt(4) = dy(t0+3*h,yy(4));\n",
+"\n",
+" yP = yy(1) + 4*h/3*(2*dydt(2)-dydt(3)+2*dydt(4));\n",
+" dydt(5) = dy(t0+4*h,yP);\n",
+" yC = yy(3) + h/3*(dydt(3)+4*dydt(4)+dydt(5));\n",
+"end\n",
+"yC = round(yC*10^4)/10^4;\n",
+"disp(yC,'y(2.0) = ')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 8.8: Milne_Predictor_Corrector_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 8.8\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f] = dy(t,y)\n",
+" f = t+y;\n",
+"endfunction\n",
+"\n",
+"\n",
+"tt = 0:0.1:0.3;\n",
+"yy = [1 1.1103 1.2428 1.3997];\n",
+"\n",
+"t0 = tt(1);\n",
+"y0 = yy(1);\n",
+"t = 2;\n",
+"h = tt(2) - tt(1);\n",
+"n = (t-t0)/h;\n",
+"for i = 1:n\n",
+" dydt(1) = dy(t0,yy(1));\n",
+" dydt(2) = dy(t0+h,yy(2));\n",
+" dydt(3) = dy(t0+2*h,yy(3));\n",
+" dydt(4) = dy(t0+3*h,yy(4));\n",
+"\n",
+" yP = yy(1) + 4*h/3*(2*dydt(2)-dydt(3)+2*dydt(4));\n",
+" dydt(5) = dy(t0+4*h,yP);\n",
+" yC = yy(3) + h/3*(dydt(3)+4*dydt(4)+dydt(5));\n",
+"end\n",
+"yC = round(yC*10^4)/10^4;\n",
+"disp(yC,'y(0.4) = ')\n",
+"\n",
+"t = [tt'; t0+4*h];\n",
+"y = [yy'; yC];\n",
+"mprintf('\n%6s %8s','t','y')\n",
+"disp([t y])"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 8.9: Adam_Moulton_Predictor_Corrector_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 8.9\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"function [f] = fun1(t,y)\n",
+" f = y - t^2;\n",
+"endfunction\n",
+"\n",
+"function [f] = rk4(t,y)\n",
+" k1 = h*fun1(t,y);\n",
+" k2 = h*fun1(t+1/2*h,y+1/2*k1);\n",
+" k3 = h*fun1(t+1/2*h,y+1/2*k2);\n",
+" k4 = h*fun1(t+h,y+k1);\n",
+" f = y + 1/6*(k1+2*k2+2*k3+k4);\n",
+"endfunction\n",
+"\n",
+"t0 = 0;\n",
+"y0 = 1;\n",
+"t = 1;\n",
+"h = 0.2;\n",
+"n = (t-t0)/h;\n",
+"y = y0;\n",
+"\n",
+"for i = 2:4\n",
+" y(i) = rk4(t0,y0);\n",
+" t0 = t0 + h;\n",
+" y0 = y(i);\n",
+"end\n",
+"\n",
+"t0 = 0;\n",
+"dydt(1) = fun1(t0,y(1));\n",
+"dydt(2) = fun1(t0+h,y(2));\n",
+"dydt(3) = fun1(t0+2*h,y(3));\n",
+"dydt(4) = fun1(t0+3*h,y(4));\n",
+"\n",
+"for i = 1:n-3\n",
+" yP = y(4) + h/24*(55*dydt(4)-59*dydt(3)+37*dydt(2)-9*dydt(1));\n",
+" dydt(5) = fun1(t0+(3+i)*h,yP);\n",
+" yC = y(4) + h/24*(9*dydt(5)+19*dydt(4)-5*dydt(3)+dydt(2));\n",
+" y = [y(2:4); yC];\n",
+" dydt = [dydt(2:4); fun1(t0+(3+i)*h,yC)]\n",
+"end\n",
+"disp(yC,'Computed Solution: y(1.0) = ')\n",
+"\n",
+"function [f] = true(t)\n",
+" f = t^2 + 2*t +2 - exp(t);\n",
+"endfunction\n",
+"ytrue = true(1.0);\n",
+"ytrue = round(ytrue*10^4)/10^4;\n",
+"disp(ytrue,'Analytical Solution: y(1.0) = ')"
+ ]
+ }
+],
+"metadata": {
+ "kernelspec": {
+ "display_name": "Scilab",
+ "language": "scilab",
+ "name": "scilab"
+ },
+ "language_info": {
+ "file_extension": ".sce",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "scilab",
+ "version": "0.7.1"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/9-Parabolic_Partial_Differential_Equations.ipynb b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/9-Parabolic_Partial_Differential_Equations.ipynb
new file mode 100644
index 0000000..173fe26
--- /dev/null
+++ b/Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao/9-Parabolic_Partial_Differential_Equations.ipynb
@@ -0,0 +1,429 @@
+{
+"cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Chapter 9: Parabolic Partial Differential Equations"
+ ]
+ },
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 9.1: Taylor_Series_Expansion.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"// Example 9.1\n",
+"// This is an analytical problem and need not be coded."
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 9.2: Initial_Boundary_Value_Problem_using_Explicit_Finite_Difference_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 9.2\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"delx = 0.1;\n",
+"delt = 0.002;\n",
+"xf = 1;\n",
+"tf = 0.006;\n",
+"x = 0:delx:xf;\n",
+"t = 0:delt:tf;\n",
+"m = length(x);\n",
+"n = length(t);\n",
+"lamda = delt/delx^2;\n",
+"\n",
+"y = zeros(n,m);\n",
+"y(1:n,1) = 0;\n",
+"y(1:n,m) = 0;\n",
+"y(1,1:m) = sin(%pi*x);\n",
+"for k = 2:n\n",
+" M1 = zeros(m-2);\n",
+" M2 = zeros(m-2,1);\n",
+" for i = 1:m-2\n",
+" M1(i,i) = 1+2*lamda;\n",
+" if i==1\n",
+" M1(i,i+1) = -lamda;\n",
+" M2(i) = y(k-1,i+1) + lamda*y(k,i);\n",
+" elseif i==m-2\n",
+" M1(i,i-1) = -lamda;\n",
+" M2(i) = y(k-1,i+1) + lamda*y(k,i+2);\n",
+" else\n",
+" M1(i,i+1) = -lamda;\n",
+" M1(i,i-1) = -lamda;\n",
+" M2(i) = y(k-1,i+1);\n",
+" end\n",
+" end\n",
+" y(k,2:m-1) = (M1\M2)';\n",
+"end\n",
+"y = round(y*10^4)/10^4;\n",
+"mprintf('%4s %7s %9s %8s %9s %9s %9s %9s %9s %9s %9s %9s %9s','n','t','x = 0.0','x = 0.1','x = 0.2','x = 0.3','x = 0.4','x = 0.5','x = 0.6','x = 0.7','x = 0.8','x = 0.9','x = 1.0');\n",
+"disp([(0:n-1)' t' y])\n",
+"\n",
+"disp('At t = 0.006:')\n",
+"disp(y(n,1:m),'Computed Solution:')\n",
+"Texact = exp(-%pi^2*tf)*sin(%pi*x);\n",
+"Texact = round(Texact*10^4)/10^4;\n",
+"disp(Texact,'Analytical Solution:')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 9.3: Initial_Boundary_Value_Problem_using_Explicit_Finite_Difference_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 9.3\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"delx = 0.1;\n",
+"delt = 0.001;\n",
+"xf = 0.5;\n",
+"tf = 0.003;\n",
+"x = 0:delx:xf;\n",
+"t = 0:delt:tf;\n",
+"m = length(x);\n",
+"n = length(t);\n",
+"r = delt/delx^2;\n",
+"\n",
+"T = zeros(m,n);\n",
+"T(1:m,1) = 0;\n",
+"delTxi = 0;\n",
+"delTxf = 1;\n",
+"\n",
+"for j = 1:n\n",
+" M1 = zeros(m,m);\n",
+" M2 = zeros(m,1);\n",
+" for i = 1:m\n",
+" if i == 1 then\n",
+" M1(i,i) = 1;\n",
+" M1(i,i+1) = -1;\n",
+" M2(i) = -delx * delTxi;\n",
+" elseif i == m then\n",
+" M1(i,i) = 1;\n",
+" M1(i,i-1) = -1;\n",
+" M2(i) = delx * delTxf;\n",
+" else\n",
+" M1(i,i) = 1;\n",
+" M2(i) = r*T(i+1,j) + (1-2*r) * T(i,j) + r*T(i-1,j);\n",
+" end\n",
+" end\n",
+" T(1:m,j+1) = (M1\M2);\n",
+"end\n",
+"T = T(:,2:n+1);\n",
+"mprintf('%4s %7s %9s %5s %7s %9s %9s %9s','n','t','x = 0.0','x=0.1','x = 0.2','x = 0.3','x = 0.4','x = 0.5');\n",
+"disp([(0:n-1)' t' T'])"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 9.4: Crank_Nicolson_Finite_Difference_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 9.4\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"delx = 0.25;\n",
+"delt = 1/32;\n",
+"xf = 1;\n",
+"tf = delt;\n",
+"x = 0:delx:xf;\n",
+"t = 0:delt:tf;\n",
+"m = length(x);\n",
+"n = length(t);\n",
+"r = delt/delx^2;\n",
+"\n",
+"\n",
+"T = zeros(m,n);\n",
+"T(1:m,1) = 1;\n",
+"T(1,1:n) = 0;\n",
+"T(m,1:n) = 0;\n",
+"\n",
+"for j = 1:n-1\n",
+" M1 = zeros(m-2,m-2);\n",
+" M2 = zeros(m-2,1);\n",
+" for i = 2:m-1\n",
+" if i == 2 then\n",
+" M1(i-1,i-1) = -2*(1+r);\n",
+" M1(i-1,i) = r;\n",
+" M2(i-1) = -(r*T(i+1,j) + 2*(1-r)*T(i,j) + r*T(i-1,j));\n",
+" elseif i == m-1\n",
+" M1(i-1,i-2) = r;\n",
+" M1(i-1,i-1) = -2*(1+r);\n",
+" M2(i-1) = -(r*T(i+1,j) + 2*(1-r)*T(i,j) + r*T(i-1,j));\n",
+" else\n",
+" M1(i-1,i-2) = r;\n",
+" M1(i-1,i-1) = -2*(1+r);\n",
+" M1(i-1,i) = r;\n",
+" M2(i-1) = -(r*T(i+1,j) + 2*(1-r)*T(i,j) + r*T(i-1,j));\n",
+" end\n",
+" end\n",
+" T(2:m-1,j+1) = (M1\M2);\n",
+"end\n",
+"T1 = M1\M2;\n",
+"for i = 1:length(T1)\n",
+" disp(strcat(['T',string(i),' = ',string(T1(i))]));\n",
+"end"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 9.5: Crank_Nicolson_Finite_Difference_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 9.5\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"delx = 1;\n",
+"delt = 1.893;\n",
+"alpha = 0.132;\n",
+"xf = 4;\n",
+"tf = delt;\n",
+"x = 0:delx:xf;\n",
+"t = 0:delt:tf;\n",
+"m = length(x);\n",
+"n = length(t);\n",
+"r = alpha*delt/delx^2;\n",
+"r = round(r*10^2)/10^2;\n",
+"\n",
+"T = zeros(m,n);\n",
+"T(1:m,1) = 1000;\n",
+"\n",
+"for j = 1:n-1\n",
+" M1 = zeros(m,m);\n",
+" M2 = zeros(m,1);\n",
+" for i = 1:m\n",
+" if i == 1 then\n",
+" M1(i,i) = -(2+2.15*r);\n",
+" M1(i,i+1) = 2*r;\n",
+" M2(i) = -(2*r*T(2,j) + (2-2.15*r)*T(1,j) + 21*r);\n",
+" elseif i == m then\n",
+" M1(i,i) = -(2+1.75*r);\n",
+" M1(i,i-1) = 2*r;\n",
+" M2(i) = -(2*r*T(m-1,j) + (2-1.75*r)*T(m,j) - 35*r);\n",
+" else\n",
+" M1(i,i-1) = r;\n",
+" M1(i,i) = -2*(1+r);\n",
+" M1(i,i+1) = r;\n",
+" M2(i) = -(r*T(i+1,j) + 2*(1-r)*T(i,j) + r*T(i-1,j));\n",
+" end\n",
+" end\n",
+" T(1:m,j+1) = (M1\M2);\n",
+"end\n",
+"disp(M1,'Coefficient Matrix:')\n",
+"disp(M2,'Constant Matrix:')\n",
+"T = round(T*10^4)/10^4;\n",
+"disp(T','Table:')"
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 9.6: Crank_Nicolson_Scheme_for_Diffusion_Equation.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"// Example 9.6\n",
+"// This is an analytical problem and need not be coded."
+ ]
+ }
+,
+{
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example 9.7: Alternate_Direction_Implicit_Method.sce"
+ ]
+ },
+ {
+"cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+"source": [
+"//Example 9.7\n",
+"\n",
+"clc\n",
+"clear\n",
+"\n",
+"h = 2;\n",
+"delt = 4;\n",
+"tf = 8;\n",
+"xf = 8;\n",
+"yf = 6;\n",
+"x = 0:h:xf;\n",
+"y = 0:h:yf;\n",
+"t = 0:delt:tf;\n",
+"m = length(x);\n",
+"n = length(y);\n",
+"p = length(t);\n",
+"r = delt/h^2;\n",
+"r = round(r*10^2)/10^2;\n",
+"\n",
+"T = 50*ones(n,m);\n",
+"T0 = T;\n",
+"T(1,1:m) = 110:-10:70;\n",
+"T(n,1:m) = 0:10:40;\n",
+"T(2:n-1,1) = [65; 25];\n",
+"T(2:n-1,m) = [60; 50];\n",
+"\n",
+"u = (m-2)*(n-2);\n",
+"index = [repmat(2:m-1,1,n-2); gsort(repmat(2:n-1,1,m-2))];\n",
+"\n",
+"M1 = zeros(u,u);\n",
+"M2 = zeros(u,1);\n",
+"for j = 2:m-1\n",
+" for i = 2:n-1\n",
+" ind = find(index(1,:)== j& index(2,:) == i);\n",
+" if j == 2 then\n",
+" M1(ind,ind) = 1+2*r;\n",
+" M1(ind,ind+1) = -r;\n",
+" M2(ind) = r*T(i,j-1) + r*T0(i-1,j) + (1-2*r)*T0(i,j) + r*T0(i+1,j);\n",
+" elseif j == m-1 then\n",
+" M1(ind,ind-1) = -r;\n",
+" M1(ind,ind) = 1+2*r;\n",
+" M2(ind) = r*T(i,j+1) + r*T0(i-1,j) + (1-2*r)*T0(i,j) + r*T0(i+1,j);\n",
+" else\n",
+" M1(ind,ind-1) = -r;\n",
+" M1(ind,ind) = 1+2*r;\n",
+" M1(ind,ind+1) = -r;\n",
+" M2(ind) = r*T0(i-1,j) + (1-2*r)*T0(i,j) + r*T0(i+1,j);\n",
+" end\n",
+" end\n",
+"end\n",
+"value = M1\M2;\n",
+"value = round(value*10^4)/10^4;\n",
+"for i = 1:length(index(1,:))\n",
+" t = index(:,i);\n",
+" T(t(2),t(1)) = value(i);\n",
+"end\n",
+"disp(T,'At t = 4:')\n",
+"T0 = T;\n",
+"\n",
+"index = gsort(index,'lc','i');\n",
+"M1 = zeros(u,u);\n",
+"M2 = zeros(u,1);\n",
+"for j = 2:m-1\n",
+" for i = 2:n-1\n",
+" ind = find(index(1,:)== j& index(2,:) == i);\n",
+" if i == 2 then\n",
+" M1(ind,ind) = 1+2*r;\n",
+" M1(ind,ind+1) = -r;\n",
+" M2(ind) = r*T(i-1,j) + r*T0(i,j-1) + (1-2*r)*T0(i,j) + r*T0(i,j+1);\n",
+" elseif i == n-1 then\n",
+" M1(ind,ind-1) = -r;\n",
+" M1(ind,ind) = 1+2*r;\n",
+" M2(ind) = r*T(i+1,j) + r*T0(i,j-1) + (1-2*r)*T0(i,j) + r*T0(i,j+1);\n",
+" else\n",
+" M1(ind,ind-1) = -r;\n",
+" M1(ind,ind) = 1+2*r;\n",
+" M1(ind,ind+1) = -r;\n",
+" M2(ind) = + r*T0(i,j-1) + (1-2*r)*T0(i,j) + r*T0(i,j+1);\n",
+" end\n",
+" end\n",
+"end\n",
+"value = M1\M2;\n",
+"value = round(value*10^4)/10^4;\n",
+"for i = 1:length(index(1,:))\n",
+" t = index(:,i);\n",
+" T(t(2),t(1)) = value(i);\n",
+"end\n",
+"disp(T,'At t = 8:')"
+ ]
+ }
+],
+"metadata": {
+ "kernelspec": {
+ "display_name": "Scilab",
+ "language": "scilab",
+ "name": "scilab"
+ },
+ "language_info": {
+ "file_extension": ".sce",
+ "help_links": [
+ {
+ "text": "MetaKernel Magics",
+ "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
+ }
+ ],
+ "mimetype": "text/x-octave",
+ "name": "scilab",
+ "version": "0.7.1"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}