diff options
author | Prashant S | 2020-04-14 10:25:32 +0530 |
---|---|---|
committer | GitHub | 2020-04-14 10:25:32 +0530 |
commit | 06b09e7d29d252fb2f5a056eeb8bd1264ff6a333 (patch) | |
tree | 2b1df110e24ff0174830d7f825f43ff1c134d1af /Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao | |
parent | abb52650288b08a680335531742a7126ad0fb846 (diff) | |
parent | 476705d693c7122d34f9b049fa79b935405c9b49 (diff) | |
download | all-scilab-tbc-books-ipynb-master.tar.gz all-scilab-tbc-books-ipynb-master.tar.bz2 all-scilab-tbc-books-ipynb-master.zip |
Initial commit
Diffstat (limited to 'Numerical_Methods_For_Scientists_And_Engineers_by_K_S_Rao')
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 +} |