From b1f5c3f8d6671b4331cef1dcebdf63b7a43a3a2b Mon Sep 17 00:00:00 2001 From: priyanka Date: Wed, 24 Jun 2015 15:03:17 +0530 Subject: initial commit / add all books --- 50/CH2/EX2.1/2_1.PNG | Bin 0 -> 13241 bytes 50/CH2/EX2.1/ex_1.sce | 18 +++++++++ 50/CH2/EX2.11/ex_11.sce | 25 +++++++++++++ 50/CH2/EX2.12/2_12.PNG | Bin 0 -> 9091 bytes 50/CH2/EX2.12/ex_12.sce | 24 ++++++++++++ 50/CH2/EX2.13/2_13.PNG | Bin 0 -> 11143 bytes 50/CH2/EX2.13/ex_13.sce | 33 +++++++++++++++++ 50/CH2/EX2.14/2_14.PNG | Bin 0 -> 8976 bytes 50/CH2/EX2.14/ex_14.sce | 30 +++++++++++++++ 50/CH2/EX2.15/2_15.PNG | Bin 0 -> 9765 bytes 50/CH2/EX2.15/ex_15.sce | 35 ++++++++++++++++++ 50/CH2/EX2.16/2_16.PNG | Bin 0 -> 11143 bytes 50/CH2/EX2.16/ex_16.sce | 38 +++++++++++++++++++ 50/CH2/EX2.17/ex_17.sce | 34 +++++++++++++++++ 50/CH2/EX2.2/2_2.PNG | Bin 0 -> 12003 bytes 50/CH2/EX2.2/ex_2.sce | 17 +++++++++ 50/CH2/EX2.23/ex_23.sce | 48 ++++++++++++++++++++++++ 50/CH2/EX2.24.1/2_24_1.PNG | Bin 0 -> 11049 bytes 50/CH2/EX2.24.1/ex_24_1.sce | 39 ++++++++++++++++++++ 50/CH2/EX2.24.2/2_24_2.PNG | Bin 0 -> 9894 bytes 50/CH2/EX2.24.2/ex_24_2.sce | 39 ++++++++++++++++++++ 50/CH2/EX2.25/2_25.PNG | Bin 0 -> 10589 bytes 50/CH2/EX2.25/ex_25.sce | 36 ++++++++++++++++++ 50/CH2/EX2.26/2_26.PNG | Bin 0 -> 9944 bytes 50/CH2/EX2.26/ex_26.sce | 36 ++++++++++++++++++ 50/CH2/EX2.27/2_27.PNG | Bin 0 -> 11063 bytes 50/CH2/EX2.27/ex_27.sce | 46 +++++++++++++++++++++++ 50/CH2/EX2.3/2_3.PNG | Bin 0 -> 11219 bytes 50/CH2/EX2.3/ex_3.sce | 35 ++++++++++++++++++ 50/CH2/EX2.4/2_4.PNG | Bin 0 -> 10551 bytes 50/CH2/EX2.4/ex_4.sce | 32 ++++++++++++++++ 50/CH2/EX2.5/ex_5.sce | 48 ++++++++++++++++++++++++ 50/CH2/EX2.6/2_6.PNG | Bin 0 -> 10673 bytes 50/CH2/EX2.6/ex_6.sce | 55 ++++++++++++++++++++++++++++ 50/CH2/EX2.7/2_7.PNG | Bin 0 -> 10955 bytes 50/CH2/EX2.7/ex_7.sce | 30 +++++++++++++++ 50/CH2/EX2.8/2_8.PNG | Bin 0 -> 10361 bytes 50/CH2/EX2.8/ex_8.sce | 31 ++++++++++++++++ 50/CH2/EX2.9/2_9.PNG | Bin 0 -> 10478 bytes 50/CH2/EX2.9/ex_9.sce | 37 +++++++++++++++++++ 50/CH3/EX3.1/ex_3_1.sce | 14 +++++++ 50/CH3/EX3.10/ex_3_10.sce | 13 +++++++ 50/CH3/EX3.11/ex_3_11.sce | 12 ++++++ 50/CH3/EX3.12/ex_3_12.sce | 28 ++++++++++++++ 50/CH3/EX3.13/ex_3_13.sce | 18 +++++++++ 50/CH3/EX3.14/ex_3_14.sce | 12 ++++++ 50/CH3/EX3.15/ex_3_15.sce | 19 ++++++++++ 50/CH3/EX3.2/ex_3_2.sce | 13 +++++++ 50/CH3/EX3.21/ex_3_21.sce | 14 +++++++ 50/CH3/EX3.22/ex_3_22.sce | 14 +++++++ 50/CH3/EX3.27/ex_3_27.sce | 24 ++++++++++++ 50/CH3/EX3.4/ex_3_4.sce | 30 +++++++++++++++ 50/CH3/EX3.5/ex_3_5.sce | 9 +++++ 50/CH3/EX3.6/ex_3_6.sce | 10 +++++ 50/CH3/EX3.8/ex_3_8.sce | 8 ++++ 50/CH3/EX3.9/ex_3_9.sce | 15 ++++++++ 50/CH4/EX4.15/ex_4_15.sce | 26 +++++++++++++ 50/CH4/EX4.20/ex_4_20.sce | 14 +++++++ 50/CH4/EX4.21/ex_4_21.sce | 17 +++++++++ 50/CH4/EX4.22/ex_4_22.sce | 30 +++++++++++++++ 50/CH4/EX4.23/ex_4_23.sce | 26 +++++++++++++ 50/CH4/EX4.3/ex_4_3.sce | 19 ++++++++++ 50/CH4/EX4.31/ex_4_31.sce | 23 ++++++++++++ 50/CH4/EX4.32/ex_4_32.sce | 45 +++++++++++++++++++++++ 50/CH4/EX4.34/ex_4_34.sce | 9 +++++ 50/CH4/EX4.35/ex_4_35.sce | 7 ++++ 50/CH4/EX4.36/ex_4_36.sce | 31 ++++++++++++++++ 50/CH4/EX4.37/ex_4_37.sce | 30 +++++++++++++++ 50/CH4/EX4.38/ex_4_38.sce | 26 +++++++++++++ 50/CH4/EX4.39/ex_4_39.sce | 34 +++++++++++++++++ 50/CH4/EX4.4/ex_4_4.sce | 12 ++++++ 50/CH4/EX4.41/ex_4_41.sce | 20 ++++++++++ 50/CH4/EX4.6/ex_4_6.sce | 18 +++++++++ 50/CH4/EX4.7/ex_4_7.sce | 15 ++++++++ 50/CH4/EX4.8/ex_4_8.sce | 15 ++++++++ 50/CH4/EX4.9/ex_4_9.sce | 16 ++++++++ 50/CH5/EX5.1/ex_5_1.sce | 12 ++++++ 50/CH5/EX5.10/ex_5_10.sce | 15 ++++++++ 50/CH5/EX5.11/ex_5_11.sce | 18 +++++++++ 50/CH5/EX5.12/ex_5_12.sce | 17 +++++++++ 50/CH5/EX5.13/ex_5_13.sce | 12 ++++++ 50/CH5/EX5.15/ex_5_15.sce | 26 +++++++++++++ 50/CH5/EX5.16/ex_5_16.sce | 16 ++++++++ 50/CH5/EX5.17/ex_5_17.sce | 24 ++++++++++++ 50/CH5/EX5.18/ex_5_18.sce | 23 ++++++++++++ 50/CH5/EX5.2/ex_5_2.sce | 9 +++++ 50/CH5/EX5.20/ex_5_20.sce | 16 ++++++++ 50/CH5/EX5.21/ex_5_21.sce | 16 ++++++++ 50/CH5/EX5.22/ex_5_22.sce | 20 ++++++++++ 50/CH5/EX5.26/ex_5_26.sce | 42 +++++++++++++++++++++ 50/CH5/EX5.27/ex_5_27.sce | 23 ++++++++++++ 50/CH5/EX5.29/ex_5_29.sce | 12 ++++++ 50/CH5/EX5.30/ex_5_30.sce | 18 +++++++++ 50/CH6/EX6.12/ex_6_12.sce | 9 +++++ 50/CH6/EX6.13/ex_6_13.sce | 11 ++++++ 50/CH6/EX6.15/ex_6_15.sce | 20 ++++++++++ 50/CH6/EX6.17/ex_6_17.sce | 17 +++++++++ 50/CH6/EX6.18/ex_6_18.sce | 8 ++++ 50/CH6/EX6.20/ex_6_20.sce | 18 +++++++++ 50/CH6/EX6.21/ex_6_21.sce | 27 ++++++++++++++ 50/CH6/EX6.25/ex_6_25.sce | 8 ++++ 50/CH6/EX6.27/ex_6_27.sce | 32 ++++++++++++++++ 50/CH6/EX6.3/ex_6_3.sce | 17 +++++++++ 50/CH6/EX6.32/ex_6_32.sce | 25 +++++++++++++ 50/CH6/EX6.4/ex_6_4.sce | 16 ++++++++ 50/CH6/EX6.9/ex_6_9.sce | 14 +++++++ 50/CH7/EX7.1/ex_7_1.sce | 27 ++++++++++++++ 50/CH7/EX7.11/ex_7_11.sce | 32 ++++++++++++++++ 50/CH7/EX7.3/ex_7_3.sce | 29 +++++++++++++++ 50/CH7/EX7.4/ex_7_4.sce | 28 ++++++++++++++ 50/CH7/EX7.5/ex_7_5.sce | 52 ++++++++++++++++++++++++++ 50/CH7/EX7.6/ex_7_6.sce | 32 ++++++++++++++++ 50/DEPENDENCIES/Euler1.sce | 26 +++++++++++++ 50/DEPENDENCIES/F28.sce | 8 ++++ 50/DEPENDENCIES/LandU.sce | 12 ++++++ 50/DEPENDENCIES/NBDP.sci | 47 ++++++++++++++++++++++++ 50/DEPENDENCIES/NDDinterpol.sci | 17 +++++++++ 50/DEPENDENCIES/NDDinterpol2.sci | 19 ++++++++++ 50/DEPENDENCIES/RK4.sci | 30 +++++++++++++++ 50/DEPENDENCIES/Vaitken.sce | 9 +++++ 50/DEPENDENCIES/Vbisection.sce | 28 ++++++++++++++ 50/DEPENDENCIES/Vbisection5.sce | 27 ++++++++++++++ 50/DEPENDENCIES/Vchebyshev.sce | 18 +++++++++ 50/DEPENDENCIES/Vgausselim.sce | 44 ++++++++++++++++++++++ 50/DEPENDENCIES/Vgaussseidel.sce | 24 ++++++++++++ 50/DEPENDENCIES/Vnewton.sce | 17 +++++++++ 50/DEPENDENCIES/Vnewton4.sce | 17 +++++++++ 50/DEPENDENCIES/Vregulafalsi.sce | 12 ++++++ 50/DEPENDENCIES/Vsecant.sce | 12 ++++++ 50/DEPENDENCIES/Vsim_eulercauchy.sce | 24 ++++++++++++ 50/DEPENDENCIES/adamsbashforth3.sci | 23 ++++++++++++ 50/DEPENDENCIES/aitken.sce | 9 +++++ 50/DEPENDENCIES/aitkeninterpol.sci | 16 ++++++++ 50/DEPENDENCIES/back.sce | 16 ++++++++ 50/DEPENDENCIES/backeuler.sci | 22 +++++++++++ 50/DEPENDENCIES/chebyshev.sce | 18 +++++++++ 50/DEPENDENCIES/cholesky.sce | 16 ++++++++ 50/DEPENDENCIES/comp_trapezoidal.sci | 34 +++++++++++++++++ 50/DEPENDENCIES/eigenvectors.sci | 37 +++++++++++++++++++ 50/DEPENDENCIES/eulermidpoint.sci | 26 +++++++++++++ 50/DEPENDENCIES/fact.sci | 6 +++ 50/DEPENDENCIES/fore.sce | 12 ++++++ 50/DEPENDENCIES/geigenvectors.sci | 37 +++++++++++++++++++ 50/DEPENDENCIES/generaliteration.sce | 21 +++++++++++ 50/DEPENDENCIES/generaliteration2.sce | 22 +++++++++++ 50/DEPENDENCIES/hermiteinterpol.sci | 37 +++++++++++++++++++ 50/DEPENDENCIES/heun.sci | 28 ++++++++++++++ 50/DEPENDENCIES/iteratedinterpol.sci | 18 +++++++++ 50/DEPENDENCIES/jacob28.sce | 9 +++++ 50/DEPENDENCIES/jacobianmat.sci | 8 ++++ 50/DEPENDENCIES/jacobiiteration.sce | 23 ++++++++++++ 50/DEPENDENCIES/jordan.sce | 20 ++++++++++ 50/DEPENDENCIES/lagrangefundamentalpoly.sci | 30 +++++++++++++++ 50/DEPENDENCIES/legrangeinterpol.sci | 18 +++++++++ 50/DEPENDENCIES/linearinterpol.sci | 3 ++ 50/DEPENDENCIES/modified_newton.sce | 17 +++++++++ 50/DEPENDENCIES/modifiedeuler.sci | 28 ++++++++++++++ 50/DEPENDENCIES/muller3.sce | 32 ++++++++++++++++ 50/DEPENDENCIES/muller5.sce | 32 ++++++++++++++++ 50/DEPENDENCIES/multipoint_iteration31.sce | 14 +++++++ 50/DEPENDENCIES/multipoint_iteration33.sce | 14 +++++++ 50/DEPENDENCIES/newton2.sce | 21 +++++++++++ 50/DEPENDENCIES/newton63.sce | 17 +++++++++ 50/DEPENDENCIES/newtonNLE.sce | 31 ++++++++++++++++ 50/DEPENDENCIES/newtonrap.sce | 18 +++++++++ 50/DEPENDENCIES/pivotgausselim.sci | 27 ++++++++++++++ 50/DEPENDENCIES/quadraticapprox.sci | 22 +++++++++++ 50/DEPENDENCIES/regulafalsi.sce | 12 ++++++ 50/DEPENDENCIES/regulafalsi4.sce | 12 ++++++ 50/DEPENDENCIES/secant4.sce | 12 ++++++ 50/DEPENDENCIES/secant64.sce | 12 ++++++ 50/DEPENDENCIES/secant65.sce | 13 +++++++ 50/DEPENDENCIES/shooting.sci | 31 ++++++++++++++++ 50/DEPENDENCIES/simRK4.sci | 27 ++++++++++++++ 50/DEPENDENCIES/sim_eulercauchy.sce | 24 ++++++++++++ 50/DEPENDENCIES/simpson.sci | 3 ++ 50/DEPENDENCIES/simpson13.sci | 34 +++++++++++++++++ 50/DEPENDENCIES/simpson38.sci | 37 +++++++++++++++++++ 50/DEPENDENCIES/straightlineapprox.sce | 16 ++++++++ 50/DEPENDENCIES/taylor.sci | 3 ++ 50/DEPENDENCIES/trapezoidal.sci | 7 ++++ 181 files changed, 3609 insertions(+) create mode 100755 50/CH2/EX2.1/2_1.PNG create mode 100755 50/CH2/EX2.1/ex_1.sce create mode 100755 50/CH2/EX2.11/ex_11.sce create mode 100755 50/CH2/EX2.12/2_12.PNG create mode 100755 50/CH2/EX2.12/ex_12.sce create mode 100755 50/CH2/EX2.13/2_13.PNG create mode 100755 50/CH2/EX2.13/ex_13.sce create mode 100755 50/CH2/EX2.14/2_14.PNG create mode 100755 50/CH2/EX2.14/ex_14.sce create mode 100755 50/CH2/EX2.15/2_15.PNG create mode 100755 50/CH2/EX2.15/ex_15.sce create mode 100755 50/CH2/EX2.16/2_16.PNG create mode 100755 50/CH2/EX2.16/ex_16.sce create mode 100755 50/CH2/EX2.17/ex_17.sce create mode 100755 50/CH2/EX2.2/2_2.PNG create mode 100755 50/CH2/EX2.2/ex_2.sce create mode 100755 50/CH2/EX2.23/ex_23.sce create mode 100755 50/CH2/EX2.24.1/2_24_1.PNG create mode 100755 50/CH2/EX2.24.1/ex_24_1.sce create mode 100755 50/CH2/EX2.24.2/2_24_2.PNG create mode 100755 50/CH2/EX2.24.2/ex_24_2.sce create mode 100755 50/CH2/EX2.25/2_25.PNG create mode 100755 50/CH2/EX2.25/ex_25.sce create mode 100755 50/CH2/EX2.26/2_26.PNG create mode 100755 50/CH2/EX2.26/ex_26.sce create mode 100755 50/CH2/EX2.27/2_27.PNG create mode 100755 50/CH2/EX2.27/ex_27.sce create mode 100755 50/CH2/EX2.3/2_3.PNG create mode 100755 50/CH2/EX2.3/ex_3.sce create mode 100755 50/CH2/EX2.4/2_4.PNG create mode 100755 50/CH2/EX2.4/ex_4.sce create mode 100755 50/CH2/EX2.5/ex_5.sce create mode 100755 50/CH2/EX2.6/2_6.PNG create mode 100755 50/CH2/EX2.6/ex_6.sce create mode 100755 50/CH2/EX2.7/2_7.PNG create mode 100755 50/CH2/EX2.7/ex_7.sce create mode 100755 50/CH2/EX2.8/2_8.PNG create mode 100755 50/CH2/EX2.8/ex_8.sce create mode 100755 50/CH2/EX2.9/2_9.PNG create mode 100755 50/CH2/EX2.9/ex_9.sce create mode 100755 50/CH3/EX3.1/ex_3_1.sce create mode 100755 50/CH3/EX3.10/ex_3_10.sce create mode 100755 50/CH3/EX3.11/ex_3_11.sce create mode 100755 50/CH3/EX3.12/ex_3_12.sce create mode 100755 50/CH3/EX3.13/ex_3_13.sce create mode 100755 50/CH3/EX3.14/ex_3_14.sce create mode 100755 50/CH3/EX3.15/ex_3_15.sce create mode 100755 50/CH3/EX3.2/ex_3_2.sce create mode 100755 50/CH3/EX3.21/ex_3_21.sce create mode 100755 50/CH3/EX3.22/ex_3_22.sce create mode 100755 50/CH3/EX3.27/ex_3_27.sce create mode 100755 50/CH3/EX3.4/ex_3_4.sce create mode 100755 50/CH3/EX3.5/ex_3_5.sce create mode 100755 50/CH3/EX3.6/ex_3_6.sce create mode 100755 50/CH3/EX3.8/ex_3_8.sce create mode 100755 50/CH3/EX3.9/ex_3_9.sce create mode 100755 50/CH4/EX4.15/ex_4_15.sce create mode 100755 50/CH4/EX4.20/ex_4_20.sce create mode 100755 50/CH4/EX4.21/ex_4_21.sce create mode 100755 50/CH4/EX4.22/ex_4_22.sce create mode 100755 50/CH4/EX4.23/ex_4_23.sce create mode 100755 50/CH4/EX4.3/ex_4_3.sce create mode 100755 50/CH4/EX4.31/ex_4_31.sce create mode 100755 50/CH4/EX4.32/ex_4_32.sce create mode 100755 50/CH4/EX4.34/ex_4_34.sce create mode 100755 50/CH4/EX4.35/ex_4_35.sce create mode 100755 50/CH4/EX4.36/ex_4_36.sce create mode 100755 50/CH4/EX4.37/ex_4_37.sce create mode 100755 50/CH4/EX4.38/ex_4_38.sce create mode 100755 50/CH4/EX4.39/ex_4_39.sce create mode 100755 50/CH4/EX4.4/ex_4_4.sce create mode 100755 50/CH4/EX4.41/ex_4_41.sce create mode 100755 50/CH4/EX4.6/ex_4_6.sce create mode 100755 50/CH4/EX4.7/ex_4_7.sce create mode 100755 50/CH4/EX4.8/ex_4_8.sce create mode 100755 50/CH4/EX4.9/ex_4_9.sce create mode 100755 50/CH5/EX5.1/ex_5_1.sce create mode 100755 50/CH5/EX5.10/ex_5_10.sce create mode 100755 50/CH5/EX5.11/ex_5_11.sce create mode 100755 50/CH5/EX5.12/ex_5_12.sce create mode 100755 50/CH5/EX5.13/ex_5_13.sce create mode 100755 50/CH5/EX5.15/ex_5_15.sce create mode 100755 50/CH5/EX5.16/ex_5_16.sce create mode 100755 50/CH5/EX5.17/ex_5_17.sce create mode 100755 50/CH5/EX5.18/ex_5_18.sce create mode 100755 50/CH5/EX5.2/ex_5_2.sce create mode 100755 50/CH5/EX5.20/ex_5_20.sce create mode 100755 50/CH5/EX5.21/ex_5_21.sce create mode 100755 50/CH5/EX5.22/ex_5_22.sce create mode 100755 50/CH5/EX5.26/ex_5_26.sce create mode 100755 50/CH5/EX5.27/ex_5_27.sce create mode 100755 50/CH5/EX5.29/ex_5_29.sce create mode 100755 50/CH5/EX5.30/ex_5_30.sce create mode 100755 50/CH6/EX6.12/ex_6_12.sce create mode 100755 50/CH6/EX6.13/ex_6_13.sce create mode 100755 50/CH6/EX6.15/ex_6_15.sce create mode 100755 50/CH6/EX6.17/ex_6_17.sce create mode 100755 50/CH6/EX6.18/ex_6_18.sce create mode 100755 50/CH6/EX6.20/ex_6_20.sce create mode 100755 50/CH6/EX6.21/ex_6_21.sce create mode 100755 50/CH6/EX6.25/ex_6_25.sce create mode 100755 50/CH6/EX6.27/ex_6_27.sce create mode 100755 50/CH6/EX6.3/ex_6_3.sce create mode 100755 50/CH6/EX6.32/ex_6_32.sce create mode 100755 50/CH6/EX6.4/ex_6_4.sce create mode 100755 50/CH6/EX6.9/ex_6_9.sce create mode 100755 50/CH7/EX7.1/ex_7_1.sce create mode 100755 50/CH7/EX7.11/ex_7_11.sce create mode 100755 50/CH7/EX7.3/ex_7_3.sce create mode 100755 50/CH7/EX7.4/ex_7_4.sce create mode 100755 50/CH7/EX7.5/ex_7_5.sce create mode 100755 50/CH7/EX7.6/ex_7_6.sce create mode 100755 50/DEPENDENCIES/Euler1.sce create mode 100755 50/DEPENDENCIES/F28.sce create mode 100755 50/DEPENDENCIES/LandU.sce create mode 100755 50/DEPENDENCIES/NBDP.sci create mode 100755 50/DEPENDENCIES/NDDinterpol.sci create mode 100755 50/DEPENDENCIES/NDDinterpol2.sci create mode 100755 50/DEPENDENCIES/RK4.sci create mode 100755 50/DEPENDENCIES/Vaitken.sce create mode 100755 50/DEPENDENCIES/Vbisection.sce create mode 100755 50/DEPENDENCIES/Vbisection5.sce create mode 100755 50/DEPENDENCIES/Vchebyshev.sce create mode 100755 50/DEPENDENCIES/Vgausselim.sce create mode 100755 50/DEPENDENCIES/Vgaussseidel.sce create mode 100755 50/DEPENDENCIES/Vnewton.sce create mode 100755 50/DEPENDENCIES/Vnewton4.sce create mode 100755 50/DEPENDENCIES/Vregulafalsi.sce create mode 100755 50/DEPENDENCIES/Vsecant.sce create mode 100755 50/DEPENDENCIES/Vsim_eulercauchy.sce create mode 100755 50/DEPENDENCIES/adamsbashforth3.sci create mode 100755 50/DEPENDENCIES/aitken.sce create mode 100755 50/DEPENDENCIES/aitkeninterpol.sci create mode 100755 50/DEPENDENCIES/back.sce create mode 100755 50/DEPENDENCIES/backeuler.sci create mode 100755 50/DEPENDENCIES/chebyshev.sce create mode 100755 50/DEPENDENCIES/cholesky.sce create mode 100755 50/DEPENDENCIES/comp_trapezoidal.sci create mode 100755 50/DEPENDENCIES/eigenvectors.sci create mode 100755 50/DEPENDENCIES/eulermidpoint.sci create mode 100755 50/DEPENDENCIES/fact.sci create mode 100755 50/DEPENDENCIES/fore.sce create mode 100755 50/DEPENDENCIES/geigenvectors.sci create mode 100755 50/DEPENDENCIES/generaliteration.sce create mode 100755 50/DEPENDENCIES/generaliteration2.sce create mode 100755 50/DEPENDENCIES/hermiteinterpol.sci create mode 100755 50/DEPENDENCIES/heun.sci create mode 100755 50/DEPENDENCIES/iteratedinterpol.sci create mode 100755 50/DEPENDENCIES/jacob28.sce create mode 100755 50/DEPENDENCIES/jacobianmat.sci create mode 100755 50/DEPENDENCIES/jacobiiteration.sce create mode 100755 50/DEPENDENCIES/jordan.sce create mode 100755 50/DEPENDENCIES/lagrangefundamentalpoly.sci create mode 100755 50/DEPENDENCIES/legrangeinterpol.sci create mode 100755 50/DEPENDENCIES/linearinterpol.sci create mode 100755 50/DEPENDENCIES/modified_newton.sce create mode 100755 50/DEPENDENCIES/modifiedeuler.sci create mode 100755 50/DEPENDENCIES/muller3.sce create mode 100755 50/DEPENDENCIES/muller5.sce create mode 100755 50/DEPENDENCIES/multipoint_iteration31.sce create mode 100755 50/DEPENDENCIES/multipoint_iteration33.sce create mode 100755 50/DEPENDENCIES/newton2.sce create mode 100755 50/DEPENDENCIES/newton63.sce create mode 100755 50/DEPENDENCIES/newtonNLE.sce create mode 100755 50/DEPENDENCIES/newtonrap.sce create mode 100755 50/DEPENDENCIES/pivotgausselim.sci create mode 100755 50/DEPENDENCIES/quadraticapprox.sci create mode 100755 50/DEPENDENCIES/regulafalsi.sce create mode 100755 50/DEPENDENCIES/regulafalsi4.sce create mode 100755 50/DEPENDENCIES/secant4.sce create mode 100755 50/DEPENDENCIES/secant64.sce create mode 100755 50/DEPENDENCIES/secant65.sce create mode 100755 50/DEPENDENCIES/shooting.sci create mode 100755 50/DEPENDENCIES/simRK4.sci create mode 100755 50/DEPENDENCIES/sim_eulercauchy.sce create mode 100755 50/DEPENDENCIES/simpson.sci create mode 100755 50/DEPENDENCIES/simpson13.sci create mode 100755 50/DEPENDENCIES/simpson38.sci create mode 100755 50/DEPENDENCIES/straightlineapprox.sce create mode 100755 50/DEPENDENCIES/taylor.sci create mode 100755 50/DEPENDENCIES/trapezoidal.sci (limited to '50') diff --git a/50/CH2/EX2.1/2_1.PNG b/50/CH2/EX2.1/2_1.PNG new file mode 100755 index 000000000..a113c4de6 Binary files /dev/null and b/50/CH2/EX2.1/2_1.PNG differ diff --git a/50/CH2/EX2.1/ex_1.sce b/50/CH2/EX2.1/ex_1.sce new file mode 100755 index 000000000..fe2fd52ce --- /dev/null +++ b/50/CH2/EX2.1/ex_1.sce @@ -0,0 +1,18 @@ + // The equation 8*x^3-12*x^2-2*x+3==0 has three real roots. + // the graph of this function can be observed here. +xset('window',0); +x=-1:.01:2.5; // defining the range of x. +deff('[y]=f(x)','y=8*x^3-12*x^2-2*x+3'); //defining the cunction +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph + +title(' y = 8*x^3-12*x^2-2*x+3') + +// from the above plot we can infre that the function has roots between +// the intervals (-1,0),(0,1),(1,2). diff --git a/50/CH2/EX2.11/ex_11.sce b/50/CH2/EX2.11/ex_11.sce new file mode 100755 index 000000000..92b538cc7 --- /dev/null +++ b/50/CH2/EX2.11/ex_11.sce @@ -0,0 +1,25 @@ + // The equation x^3-5*x+1==0 has real roots. + // the graph of this function can be observed here. +xset('window',10); +x=-2:.01:4; // defining the range of x. +deff('[y]=f(x)','y=x^3-5*x+1'); //defining the cunction. +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = x^3-5*x+1') + +// from the above plot we can infre that the function has roots between +// the intervals (0,1),(2,3). +// since we have been asked for the smallest positive root of the equation, +// we are intrested on the interval (0,1) + + + //sollution by muller method to 3 iterations. + +muller3(0,.5,1,f) + diff --git a/50/CH2/EX2.12/2_12.PNG b/50/CH2/EX2.12/2_12.PNG new file mode 100755 index 000000000..4c46e2515 Binary files /dev/null and b/50/CH2/EX2.12/2_12.PNG differ diff --git a/50/CH2/EX2.12/ex_12.sce b/50/CH2/EX2.12/ex_12.sce new file mode 100755 index 000000000..42225d058 --- /dev/null +++ b/50/CH2/EX2.12/ex_12.sce @@ -0,0 +1,24 @@ + // The equation cos(x)-x*%e^x==0 has real roots. + // the graph of this function can be observed here. +xset('window',8); +x=-1:.001:2; // defining the range of x. +deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the cunction. +deff('[y]=fp(x)','y=-sin(x)-x*%e^x-%e^x'); +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = cos(x)-x*%e^x') + +// from the above plot we can infre that the function has root between +// the interval (0,1) + + + //sollution by muller method to 5 iterations. + + +muller5(-1,0,1,f) diff --git a/50/CH2/EX2.13/2_13.PNG b/50/CH2/EX2.13/2_13.PNG new file mode 100755 index 000000000..0dca7c750 Binary files /dev/null and b/50/CH2/EX2.13/2_13.PNG differ diff --git a/50/CH2/EX2.13/ex_13.sce b/50/CH2/EX2.13/ex_13.sce new file mode 100755 index 000000000..1a87dd5e8 --- /dev/null +++ b/50/CH2/EX2.13/ex_13.sce @@ -0,0 +1,33 @@ + // The equation x^3-5*x+1==0 has real roots. + // the graph of this function can be observed here. +xset('window',12); +x=-2:.01:4; // defining the range of x. +deff('[y]=f(x)','y=x^3-5*x+1'); //defining the function. +deff('[y]=fp(x)','y=3*x^2-5'); +deff('[y]=fpp(x)','y=6*x'); +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = x^3-5*x+1') + +// from the above plot we can infre that the function has roots between +// the intervals (0,1),(2,3). +// since we have been asked for the smallest positive root of the equation, +// we are intrested on the interval (0,1) +// a=0;b=1, + + +// solution by chebyshev method + +// the approximate root after 4 iterations can be observed. + + +chebyshev(0.5,f,fp) + + +// hence the approximate root witin the permissible error of 10^-15 is .2016402, \ No newline at end of file diff --git a/50/CH2/EX2.14/2_14.PNG b/50/CH2/EX2.14/2_14.PNG new file mode 100755 index 000000000..884daebb0 Binary files /dev/null and b/50/CH2/EX2.14/2_14.PNG differ diff --git a/50/CH2/EX2.14/ex_14.sce b/50/CH2/EX2.14/ex_14.sce new file mode 100755 index 000000000..54e521762 --- /dev/null +++ b/50/CH2/EX2.14/ex_14.sce @@ -0,0 +1,30 @@ + + + // The equation 1/x-7==0 has a real root. + // the graph of this function can be observed here. +xset('window',13); +x=0.001:.001:.25; // defining the range of x. +deff('[y]=f(x)','y=1/x-7'); //defining the function. +deff('[y]=fp(x)','y=-1/x^2'); +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y =1/x-7') + +// from the above plot we can infre that the function has roots between +// the interval (0,2/7) + + + //solution by chebyshev method + + + chebyshev(0.1,f,fp) //calling the pre-defined function 'chebyshev' to find the approximate root in the range of (0,2/7). + + + + \ No newline at end of file diff --git a/50/CH2/EX2.15/2_15.PNG b/50/CH2/EX2.15/2_15.PNG new file mode 100755 index 000000000..5c4de8ecc Binary files /dev/null and b/50/CH2/EX2.15/2_15.PNG differ diff --git a/50/CH2/EX2.15/ex_15.sce b/50/CH2/EX2.15/ex_15.sce new file mode 100755 index 000000000..b8bfdab9a --- /dev/null +++ b/50/CH2/EX2.15/ex_15.sce @@ -0,0 +1,35 @@ + // The equation cos(x)-x*%e^x==0 has real roots. + // the graph of this function can be observed here. +xset('window',8); +x=-1:.001:2; // defining the range of x. +deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the cunction. +deff('[y]=fp(x)','y=-sin(x)-x*%e^x-%e^x'); +deff('[y]=fpp(x)','y=-cos(x)-x*%e^x-2*%e^x'); +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = cos(x)-x*%e^x') + +// from the above plot we can infre that the function has root between +// the interval (0,1) + + +// a=0;b=1, + + + + // solution by chebyshev with a permissible error of 10^-15. + +// we call a user-defined function 'chebyshev' so as to find the approximate +// root of the equation within the defined permissible error limit. + +chebyshev(1,f,fp) + + + +// hence the approximate root witin the permissible error of 10^-15 is \ No newline at end of file diff --git a/50/CH2/EX2.16/2_16.PNG b/50/CH2/EX2.16/2_16.PNG new file mode 100755 index 000000000..dd5cb960c Binary files /dev/null and b/50/CH2/EX2.16/2_16.PNG differ diff --git a/50/CH2/EX2.16/ex_16.sce b/50/CH2/EX2.16/ex_16.sce new file mode 100755 index 000000000..d1a65a8ca --- /dev/null +++ b/50/CH2/EX2.16/ex_16.sce @@ -0,0 +1,38 @@ + // The equation x^3-5*x+1==0 has real roots. + // the graph of this function can be observed here. +xset('window',15); +x=-2:.01:4; // defining the range of x. +deff('[y]=f(x)','y=x^3-5*x+1'); //defining the function. +deff('[y]=fp(x)','y=3*x^2-5'); +deff('[y]=fpp(x)','y=6*x'); +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = x^3-5*x+1') + +// from the above plot we can infre that the function has roots between +// the intervals (0,1),(2,3). +// since we have been asked for the smallest positive root of the equation, +// we are intrested on the interval (0,1) +// a=0;b=1, + + +// solution by multipoint iteration method + +// the approximate root after 3 iterations can be observed. + + +multipoint_iteration31(0.5,f,fp) + +// hence the approximate root witin the permissible error of 10^-15 is .201640, + + + +multipoint_iteration33(0.5,f,fp) + +// hence the approximate root witin the permissible error of 10^-15 is .201640, \ No newline at end of file diff --git a/50/CH2/EX2.17/ex_17.sce b/50/CH2/EX2.17/ex_17.sce new file mode 100755 index 000000000..f818a0af8 --- /dev/null +++ b/50/CH2/EX2.17/ex_17.sce @@ -0,0 +1,34 @@ + // The equation cos(x)-x*%e^x==0 has real roots. + // the graph of this function can be observed here. +xset('window',8); +x=-1:.001:2; // defining the range of x. +deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the function. +deff('[y]=fp(x)','y=-sin(x)-x*%e^x-%e^x'); +deff('[y]=fpp(x)','y=-cos(x)-x*%e^x-2*%e^x'); +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = cos(x)-x*%e^x') + +// from the above plot we can infre that the function has root between +// the interval (0,1) + + +// a=0;b=1, + + + + // solution by multipoint_iteration method using the formula given in equation no.2.33. + +// we call a user-defined function 'multipoint_iteration33' so as to find the approximate +// root of the equation within the defined permissible error limit. + +multipoint_iteration33(1,f,fp) + + +// hence the approximate root witin the permissible error of 10^-5 is 0.5177574. \ No newline at end of file diff --git a/50/CH2/EX2.2/2_2.PNG b/50/CH2/EX2.2/2_2.PNG new file mode 100755 index 000000000..c38387579 Binary files /dev/null and b/50/CH2/EX2.2/2_2.PNG differ diff --git a/50/CH2/EX2.2/ex_2.sce b/50/CH2/EX2.2/ex_2.sce new file mode 100755 index 000000000..3fa2a9a1b --- /dev/null +++ b/50/CH2/EX2.2/ex_2.sce @@ -0,0 +1,17 @@ + // The equation cos(x)-x*%e^x==0 has real roots. + // the graph of this function can be observed here. +xset('window',1); +x=0:.01:2; // defining the range of x. +deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the cunction. +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = cos(x)-x*%e^x') + +// from the above plot we can infre that the function has root between +// the interval (0,1) diff --git a/50/CH2/EX2.23/ex_23.sce b/50/CH2/EX2.23/ex_23.sce new file mode 100755 index 000000000..3a72496d0 --- /dev/null +++ b/50/CH2/EX2.23/ex_23.sce @@ -0,0 +1,48 @@ + // The equation 3*x^3+4*x^2+4*x+1==0 has three real roots. + // the graph of this function can be observed here. +xset('window',22); +x=-1.5:.001:1.5; // defining the range of x. +deff('[y]=f(x)','y=3*x^3+4*x^2+4*x+1'); //defining the cunction +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph + +title(' y =3*x^3+4*x^2+4*x+1') + +// from the above plot we can infre that the function has root between +// the interval (-1,0), + +x0=-.5; // initial approximation + + +// let the iterative function g(x) be x+A*(3*x^3+4*x^2+4*x+1) =g(x); + +// gp(x)=(1+A*(9*x^2+8*x+4) ) +// we need to choose a value for A , which makes abs(gp(x0))<1 + +// hence abs(gp(x0))=abs(1+9*A/4) + +A=-1:.1:1; + +abs(1+9*A/4) // tryin to check the values of abs(gp(x0)) for different values of A. + + +// from the above values of 'A' and the values of 'abs(gp(x0))', +// we can infer that for the vales of 'A 'in the range (-.8,0) g(x ) will be giving a converging solution , + +// hence deliberatele we choose a to be -0.5, + +A=-0.5; + +deff('[y]=g(x)','y= x-0.5*(3*x^3+4*x^2+4*x+1)'); +deff('[y]=gp(x)','y= 1-0.5*(9*x^2+8*x+4)'); // hence defining g(x) and gp(x), +generaliteration(x0,g,gp) + + + + diff --git a/50/CH2/EX2.24.1/2_24_1.PNG b/50/CH2/EX2.24.1/2_24_1.PNG new file mode 100755 index 000000000..8a9e4a241 Binary files /dev/null and b/50/CH2/EX2.24.1/2_24_1.PNG differ diff --git a/50/CH2/EX2.24.1/ex_24_1.sce b/50/CH2/EX2.24.1/ex_24_1.sce new file mode 100755 index 000000000..5745698d7 --- /dev/null +++ b/50/CH2/EX2.24.1/ex_24_1.sce @@ -0,0 +1,39 @@ + // The equation x^3-5*x+1==0 has real roots. + // the graph of this function can be observed here. +xset('window',2); +x=-2:.01:4; // defining the range of x. +deff('[y]=f(x)','y=x^3-5*x+1'); //defining the function. +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = x^3-5*x+1') + +// from the above plot we can infre that the function has roots between +// the intervals (0,1),(2,3). +// since we have been asked for the smallest positive root of the equation, +// we are intrested on the interval (0,1) + +x0=.5; + + //solution using linear iteration method for the first two iterations and aitken's process two times for the third iteration. + + deff('[y]=g(x)','y=1/5*(x^3+1)'); + deff('[y]=gp(x)','y=1/5*(3*x^2)'); + + +generaliteration2(x0,g,gp) + + +// from the above iterations performed we can infer that- +x1=0.225; +x2=0.202278; + + + + +aitken(x0,x1,x2,g) // calling the aitken method for one iteration \ No newline at end of file diff --git a/50/CH2/EX2.24.2/2_24_2.PNG b/50/CH2/EX2.24.2/2_24_2.PNG new file mode 100755 index 000000000..1b89face5 Binary files /dev/null and b/50/CH2/EX2.24.2/2_24_2.PNG differ diff --git a/50/CH2/EX2.24.2/ex_24_2.sce b/50/CH2/EX2.24.2/ex_24_2.sce new file mode 100755 index 000000000..9ff6fe046 --- /dev/null +++ b/50/CH2/EX2.24.2/ex_24_2.sce @@ -0,0 +1,39 @@ + // The equation x-%e^-x==0 has real roots. + // the graph of this function can be observed here. +xset('window',24); +x=-3:.01:4; // defining the range of x. +deff('[y]=f(x)','y=x-%e^-x'); //defining the cunction. +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = x-%e^-x') + +// from the above plot we can infre that the function has root between +// the interval (0,1) + +x0=1; + + //solution using linear iteration method for the first two iterations and aitken's process two times for the third iteration. + + + + deff('[y]=g(x)','y=%e^-x'); + deff('[y]=gp(x)','y=-%e^-x'); + + +generaliteration2(x0,g,gp) + + +// from the above iterations performed we can infer that- +x1=0.367879; +x2=0.692201; + + + + +aitken(x0,x1,x2,g) // calling the aitken method for one iteration \ No newline at end of file diff --git a/50/CH2/EX2.25/2_25.PNG b/50/CH2/EX2.25/2_25.PNG new file mode 100755 index 000000000..1478a1d27 Binary files /dev/null and b/50/CH2/EX2.25/2_25.PNG differ diff --git a/50/CH2/EX2.25/ex_25.sce b/50/CH2/EX2.25/ex_25.sce new file mode 100755 index 000000000..a9c110187 --- /dev/null +++ b/50/CH2/EX2.25/ex_25.sce @@ -0,0 +1,36 @@ + // The equation cos(x)-x*%e^x==0 has real roots. + // the graph of this function can be observed here. +xset('window',25); +x=0:.01:2; // defining the range of x. +deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the cunction. +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = cos(x)-x*%e^x') + +// from the above plot we can infre that the function has root between +// the interval (0,1) + +x0=0; + + //solution using linear iteration method for the first two iterations and aitken's process two times for the third iteration. + + deff('[y]=g(x)','y=x+1/2*(cos(x)-x*%e^x)'); + deff('[y]=gp(x)','y=1+1/2*(-sin(x)-x*%e^x-%e^x)'); + + +generaliteration2(x0,g,gp) + + +// from the above iterations performed we can infer that- +x1=0.50000000; +x2=0.5266110; + + + +aitken(x0,x1,x2,g) // calling the aitken method for one iteration diff --git a/50/CH2/EX2.26/2_26.PNG b/50/CH2/EX2.26/2_26.PNG new file mode 100755 index 000000000..d303cd438 Binary files /dev/null and b/50/CH2/EX2.26/2_26.PNG differ diff --git a/50/CH2/EX2.26/ex_26.sce b/50/CH2/EX2.26/ex_26.sce new file mode 100755 index 000000000..2e5761bf3 --- /dev/null +++ b/50/CH2/EX2.26/ex_26.sce @@ -0,0 +1,36 @@ + // The equation x^3-7*x^2+16*x-12==0 has real roots. + // the graph of this function can be observed here. +xset('window',25); +x=0:.001:4; // defining the range of x. +deff('[y]=f(x)','y=x^3-7*x^2+16*x-12'); //defining the cunction. +deff('[y]=fp(x)','y=3*x^2-14*x+16'); +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = x^3-7*x^2+16*x-12') + + + + +// given that the equation has double roots at x=2 hence m=2; + +m=2; + + // solution by newton raphson method + + +newton(1,f,fp) // calling the user defined function + + + + + //solution by modified newton raphsons mathod + + + +modified_newton(1,f,fp) \ No newline at end of file diff --git a/50/CH2/EX2.27/2_27.PNG b/50/CH2/EX2.27/2_27.PNG new file mode 100755 index 000000000..8cef9c44f Binary files /dev/null and b/50/CH2/EX2.27/2_27.PNG differ diff --git a/50/CH2/EX2.27/ex_27.sce b/50/CH2/EX2.27/ex_27.sce new file mode 100755 index 000000000..2d97164c5 --- /dev/null +++ b/50/CH2/EX2.27/ex_27.sce @@ -0,0 +1,46 @@ + // The equation 27*x^5+27*x^4+36*x^3+28*x^2+9*x+1==0 has real roots. + // the graph of this function can be observed here. +xset('window',26); +x=-2:.001:3; // defining the range of x. +deff('[y]=f(x)','y=27*x^5+27*x^4+36*x^3+28*x^2+9*x+1'); //defining the cunction. +deff('[y]=fp(x)','y=27*5*x^4+27*4*x^3+36*3*x^2+28*2*x+9'); +deff('[y]=fpp(x)','y=27*5*4*x^3+27*4*3*x^2+36*3*2*x+28*2'); +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = 27*x^5+27*x^4+36*x^3+28*x^2+9*x+1') + + + + + // solution by newton raphson method as per the equation no. 2.14 + + +newton(-1,f,fp) // calling the user defined function + + +newton4(-1,f,fp) + + + // solution by newton raphson method as per the equation no. 2.63 + +newton63(-1,f,fp,fpp) // calling the user defined function + // solution by the secant method defined to satisfy the equation no.2.64. + + +secant64(0,-1,f,fp) + + + + + + + // solution by the secant method defined to satisfy the equation no.2.65. + + +secant65(0,-.5,f) diff --git a/50/CH2/EX2.3/2_3.PNG b/50/CH2/EX2.3/2_3.PNG new file mode 100755 index 000000000..1d812544f Binary files /dev/null and b/50/CH2/EX2.3/2_3.PNG differ diff --git a/50/CH2/EX2.3/ex_3.sce b/50/CH2/EX2.3/ex_3.sce new file mode 100755 index 000000000..7add08896 --- /dev/null +++ b/50/CH2/EX2.3/ex_3.sce @@ -0,0 +1,35 @@ + // The equation x^3-5*x+1==0 has real roots. + // the graph of this function can be observed here. +xset('window',2); +x=-2:.01:4; // defining the range of x. +deff('[y]=f(x)','y=x^3-5*x+1'); //defining the cunction. +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = x^3-5*x+1') + +// from the above plot we can infre that the function has roots between +// the intervals (0,1),(2,3). +// since we have been asked for the smallest positive root of the equation, +// we are intrested on the interval (0,1) +// a=0;b=1, + +// we call a user-defined function 'bisection' so as to find the approximate +// root of the equation with a defined permissible error. + +bisection(0,1,f) + +// since in the example 2.3 we have been asked to perform 5 itterations +// the approximate root after 5 iterations can be observed. + + + +bisection5(0,1,f) + + +// hence the approximate root after 5 iterations is 0.203125 witin the permissible error of 10^-4, \ No newline at end of file diff --git a/50/CH2/EX2.4/2_4.PNG b/50/CH2/EX2.4/2_4.PNG new file mode 100755 index 000000000..e111288e6 Binary files /dev/null and b/50/CH2/EX2.4/2_4.PNG differ diff --git a/50/CH2/EX2.4/ex_4.sce b/50/CH2/EX2.4/ex_4.sce new file mode 100755 index 000000000..751cc0f5c --- /dev/null +++ b/50/CH2/EX2.4/ex_4.sce @@ -0,0 +1,32 @@ + // The equation cos(x)-x*%e^x==0 has real roots. + // the graph of this function can be observed here. +xset('window',3); +x=0:.01:2; // defining the range of x. +deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the cunction. +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = cos(x)-x*%e^x') + +// from the above plot we can infre that the function has root between +// the interval (0,1) + + +// a=0;b=1, + +// we call a user-defined function 'bisection' so as to find the approximate +// root of the equation with a defined permissible error. + +bisection(0,1,f) + +// since in the example 2.4 we have been asked to perform 5 itterations , + +bisection5(0,1,f) + + +// hence the approximate root after 5 iterations is 0.515625 witin the permissible error of 10^-4, \ No newline at end of file diff --git a/50/CH2/EX2.5/ex_5.sce b/50/CH2/EX2.5/ex_5.sce new file mode 100755 index 000000000..0a6ff1b8e --- /dev/null +++ b/50/CH2/EX2.5/ex_5.sce @@ -0,0 +1,48 @@ + // The equation x^3-5*x+1==0 has real roots. + // the graph of this function can be observed here. +xset('window',4); +x=-2:.01:4; // defining the range of x. +deff('[y]=f(x)','y=x^3-5*x+1'); //defining the cunction. +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = x^3-5*x+1') + +// from the above plot we can infre that the function has roots between +// the intervals (0,1),(2,3). +// since we have been given the interval to be considered as (0,1) +// a=0;b=1, + + + // Solution by secant method + + + + + +// since in the example 2.5 we have been asked to perform 4 itterations , +secant4(0,1,f) // we call a user-defined function 'bisection' so as to find the approximate +// root of the equation with a defined permissible error. + + + +// hence the approximate root occured in secant method after 4 iterations is 0.201640 witin the permissible error of 10^-4, + + + + // solution by regular falsi method + + +// since in the example 2.5 we have been asked to perform 4 itterations , + +regulafalsi4(0,1,f) // we call a user-defined function 'regularfalsi4' so as to find the approximate +// root of the equation with a defined permissible error. + + + +// hence the approximate root occured in regularfalsi method after 4 iterations is 0.201640 witin the permissible error of 10^-4, \ No newline at end of file diff --git a/50/CH2/EX2.6/2_6.PNG b/50/CH2/EX2.6/2_6.PNG new file mode 100755 index 000000000..3fafbbe75 Binary files /dev/null and b/50/CH2/EX2.6/2_6.PNG differ diff --git a/50/CH2/EX2.6/ex_6.sce b/50/CH2/EX2.6/ex_6.sce new file mode 100755 index 000000000..0aa95ce3a --- /dev/null +++ b/50/CH2/EX2.6/ex_6.sce @@ -0,0 +1,55 @@ + // The equation cos(x)-x*%e^x==0 has real roots. + // the graph of this function can be observed here. +xset('window',3); +x=0:.01:2; // defining the range of x. +deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the cunction. +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = cos(x)-x*%e^x') + +// from the above plot we can infre that the function has root between +// the interval (0,1) + + +// a=0;b=1, + + + // Solution by secant method + + + + + +// since in the example 2.6 we have no specification of the no. of itterations , +// we define a function 'secant' and execute it. + + + +secant(0,1,f) // we call a user-defined function 'secant' so as to find the approximate +// root of the equation with a defined permissible error. + + + +// hence the approximate root occured in secant method witin the permissible error of 10^-5 is , + + + + // solution by regular falsi method + + + + + +// since in the example 2.6 we have no specification of the no. of itterations , + + +regulafalsi(0,1,f) // we call a user-defined function 'regularfalsi' so as to find the approximate +// root of the equation with a defined permissible error. + + diff --git a/50/CH2/EX2.7/2_7.PNG b/50/CH2/EX2.7/2_7.PNG new file mode 100755 index 000000000..dcb0e7198 Binary files /dev/null and b/50/CH2/EX2.7/2_7.PNG differ diff --git a/50/CH2/EX2.7/ex_7.sce b/50/CH2/EX2.7/ex_7.sce new file mode 100755 index 000000000..ad7b083d7 --- /dev/null +++ b/50/CH2/EX2.7/ex_7.sce @@ -0,0 +1,30 @@ + // The equation x^3-5*x+1==0 has real roots. + // the graph of this function can be observed here. +xset('window',6); +x=-2:.01:4; // defining the range of x. +deff('[y]=f(x)','y=x^3-5*x+1'); //defining the function. +deff('[y]=fp(x)','y=3*x^2-5'); +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = x^3-5*x+1') + +// from the above plot we can infre that the function has roots between +// the intervals (0,1),(2,3). +// since we have been asked for the smallest positive root of the equation, +// we are intrested on the interval (0,1) +// a=0;b=1, + +// since in the example 2.7 we have been asked to perform 4 itterations , +// the approximate root after 4 iterations can be observed. + + +newton4(0.5,f,fp) + + +// hence the approximate root after 4 iterations is 0.201640 witin the permissible error of 10^-15, \ No newline at end of file diff --git a/50/CH2/EX2.8/2_8.PNG b/50/CH2/EX2.8/2_8.PNG new file mode 100755 index 000000000..bb4e1f32b Binary files /dev/null and b/50/CH2/EX2.8/2_8.PNG differ diff --git a/50/CH2/EX2.8/ex_8.sce b/50/CH2/EX2.8/ex_8.sce new file mode 100755 index 000000000..44ac39cef --- /dev/null +++ b/50/CH2/EX2.8/ex_8.sce @@ -0,0 +1,31 @@ + // The equation x^3-17==0 has three real roots. + // the graph of this function can be observed here. +xset('window',7); +x=-5:.001:5; // defining the range of x. +deff('[y]=f(x)','y=x^3-17'); //defining the cunction. +deff('[y]=fp(x)','y=3*x^2'); +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = x^3-17') + +// from the above plot we can infre that the function has root between +// the interval (2,3). + + + + //solution by newton raphson's method + + + + // since in example no.2.8 we have been asked to perform 4 iterations ,we define a fuction newton4'' which does newton raphson's method of finding approximate root upto 4 iterations, + + + + newton4(2,f,fp) //calling the pre-defined function 'newton4'. + \ No newline at end of file diff --git a/50/CH2/EX2.9/2_9.PNG b/50/CH2/EX2.9/2_9.PNG new file mode 100755 index 000000000..a172ae0c6 Binary files /dev/null and b/50/CH2/EX2.9/2_9.PNG differ diff --git a/50/CH2/EX2.9/ex_9.sce b/50/CH2/EX2.9/ex_9.sce new file mode 100755 index 000000000..5fdd426b7 --- /dev/null +++ b/50/CH2/EX2.9/ex_9.sce @@ -0,0 +1,37 @@ + // The equation cos(x)-x*%e^x==0 has real roots. + // the graph of this function can be observed here. +xset('window',8); +x=-1:.001:2; // defining the range of x. +deff('[y]=f(x)','y=cos(x)-x*%e^x'); //defining the cunction. +deff('[y]=fp(x)','y=-sin(x)-x*%e^x-%e^x'); +y=feval(x,f); + +a=gca(); + +a.y_location = "origin"; + +a.x_location = "origin"; +plot(x,y) // instruction to plot the graph +title(' y = cos(x)-x*%e^x') + +// from the above plot we can infre that the function has root between +// the interval (0,1) + + +// a=0;b=1, + + + + // solution by newton raphson's method with a permissible error of 10^-8. + + +// we call a user-defined function 'newton' so as to find the approximate +// root of the equation within the defined permissible error limit. + +newton(1,f,fp) + + + + + +// hence the approximate root witin the permissible error of 10^-8 is 0.5177574. \ No newline at end of file diff --git a/50/CH3/EX3.1/ex_3_1.sce b/50/CH3/EX3.1/ex_3_1.sce new file mode 100755 index 000000000..473e588b6 --- /dev/null +++ b/50/CH3/EX3.1/ex_3_1.sce @@ -0,0 +1,14 @@ +// example 3.1 +//Positive definite + +A=[12 4 -1;4 7 1;-1 1 6]; +//check if the determinant of the leading minors is a positive value- + +A1=A(1); +det(A1) +A2=A(1:2,1:2); +det(A2) +A3=A(1:3,1:3); +det(A3) + +//we observe that the determinant of the leading minors is a positive ,hence the given matrix A is a positive definite. \ No newline at end of file diff --git a/50/CH3/EX3.10/ex_3_10.sce b/50/CH3/EX3.10/ex_3_10.sce new file mode 100755 index 000000000..05f721b06 --- /dev/null +++ b/50/CH3/EX3.10/ex_3_10.sce @@ -0,0 +1,13 @@ +//example no. 3.10 +//solve system by decomposition method + +A=[1 1 1;4 3 -1;3 5 3] +n=3; + +b=[1;6;4] + +[U,L]=LandU(A,3) + +Z=fore(L,b) + +X=back(U,Z) \ No newline at end of file diff --git a/50/CH3/EX3.11/ex_3_11.sce b/50/CH3/EX3.11/ex_3_11.sce new file mode 100755 index 000000000..fdf0b3b7c --- /dev/null +++ b/50/CH3/EX3.11/ex_3_11.sce @@ -0,0 +1,12 @@ +//example no.3.11 +//caption: Inverse using LU decomposition + +A=[3 2 1;2 3 2;1 2 2] + +[U,L]=LandU(A,3) // call LandU function to evaluate U,L of A, + +//since A=L*U , +// inv(A)=inv(U)*inv(L) +// let inv(A)=AI + +AI=U^-1*L^-1 \ No newline at end of file diff --git a/50/CH3/EX3.12/ex_3_12.sce b/50/CH3/EX3.12/ex_3_12.sce new file mode 100755 index 000000000..4599b1c23 --- /dev/null +++ b/50/CH3/EX3.12/ex_3_12.sce @@ -0,0 +1,28 @@ +//example no. 3.12 +//solve system by decomposition method + +A=[1 1 -1;2 2 5;3 2 -3] + +b=[2;-3;6] + + + + // hence we can observe that LU decomposition method fails to solve this system since the pivot L(2,2)=0; + + + //we note that the coefficient matrix is not a positive definite matrix and hence its LU decomposition is not guaranteed, + + + //if we interchange the rows of A as shown below the LU decomposition would work, + + A=[3 2 -3;2 2 5;1 1 -1] + + b=[6;-3;2] + + [U,L]=LandU(A,3) // call LandU function to evaluate U,L of A, + +n=3; +Z=fore(L,b); + +X=back(U,Z) + \ No newline at end of file diff --git a/50/CH3/EX3.13/ex_3_13.sce b/50/CH3/EX3.13/ex_3_13.sce new file mode 100755 index 000000000..319071f40 --- /dev/null +++ b/50/CH3/EX3.13/ex_3_13.sce @@ -0,0 +1,18 @@ +//example no. 3.13 +//solve system by LU decomposition method + +A=[2 1 1 -2;4 0 2 1;3 2 2 0;1 3 2 -1] + +b=[-10;8;7;-5] + +[U,L]=LandU(A,4) +n=4; +Z=fore(L,b); + +X=back(U,Z) + +//since A=L*U , +// inv(A)=inv(U)*inv(L) +// let inv(A)=AI + +AI=U^-1*L^-1 diff --git a/50/CH3/EX3.14/ex_3_14.sce b/50/CH3/EX3.14/ex_3_14.sce new file mode 100755 index 000000000..3590d9c29 --- /dev/null +++ b/50/CH3/EX3.14/ex_3_14.sce @@ -0,0 +1,12 @@ +//example no. 3.14 +//solve system by cholesky method + +A=[1 2 3;2 8 22;3 22 82] + +b=[5;6;-10] + +L=cholesky (A,3) //call cholesky function to evaluate the root of the system +n=3; +Z=fore(L,b); + +X=back(L',Z) \ No newline at end of file diff --git a/50/CH3/EX3.15/ex_3_15.sce b/50/CH3/EX3.15/ex_3_15.sce new file mode 100755 index 000000000..902ce414b --- /dev/null +++ b/50/CH3/EX3.15/ex_3_15.sce @@ -0,0 +1,19 @@ +//example no. 3.15 +//solve system by cholesky method + +A=[4 -1 0 0;-1 4 -1 0;0 -1 4 -1;0 0 -1 4] + +b=[1;0;0;0] + +L=cholesky (A,4) //call cholesky function to evaluate the root of the system + +n=4; +Z=fore(L,b); + +X=back(L',Z) + +//since A=L*L' , +// inv(A)=inv(L')*inv(L) +// let inv(A)=AI + +AI=L'^-1*L^-1 \ No newline at end of file diff --git a/50/CH3/EX3.2/ex_3_2.sce b/50/CH3/EX3.2/ex_3_2.sce new file mode 100755 index 000000000..bc98fc202 --- /dev/null +++ b/50/CH3/EX3.2/ex_3_2.sce @@ -0,0 +1,13 @@ +//example no.3.2 +// show if a given matrix 'A' possesses 'property A' + +A=[2 -1 0 -1;-1 2 -1 0;0 -1 2 0;0 0 -1 2 ] + +P=[0 0 0 1;0 1 0 0;0 0 1 0;1 0 0 0] // let us take the pemutation matrix as P + +// then we find that + +P*A*P' + +// is of the form (3.2) . hence the matrix 'A' has property A + diff --git a/50/CH3/EX3.21/ex_3_21.sce b/50/CH3/EX3.21/ex_3_21.sce new file mode 100755 index 000000000..890f06c4d --- /dev/null +++ b/50/CH3/EX3.21/ex_3_21.sce @@ -0,0 +1,14 @@ +//example no. 3.21 +//solve the system by jacobi iteration method + +A=[4 1 1;1 5 2;1 2 3] + +b=[2 ;-6;-4] + +N=3; //no. of ierations +n=3; // order of the matrix is n*n + +X=[.5;-.5;-.5] //initial approximation + + +jacobiiteration(A,n,N,X,b) //call the function which performs jacobi iteration method to solve the system \ No newline at end of file diff --git a/50/CH3/EX3.22/ex_3_22.sce b/50/CH3/EX3.22/ex_3_22.sce new file mode 100755 index 000000000..4f6e78114 --- /dev/null +++ b/50/CH3/EX3.22/ex_3_22.sce @@ -0,0 +1,14 @@ +//example no. 3.22 +//solve the system by gauss seidel method + +A=[2 -1 0;-1 2 -1;0 -1 2] + +b=[7;1;1] + +N=3; //no. of ierations +n=3; // order of the matrix is n*n + +X=[0;0;0] //initial approximation + + +gaussseidel(A,n,N,X,b) //call the function which performs gauss seidel method to solve the system \ No newline at end of file diff --git a/50/CH3/EX3.27/ex_3_27.sce b/50/CH3/EX3.27/ex_3_27.sce new file mode 100755 index 000000000..242a16db7 --- /dev/null +++ b/50/CH3/EX3.27/ex_3_27.sce @@ -0,0 +1,24 @@ +// example 3.27 +// a) find eigenvalue and eigen vector; +// b) verify inv(S)*A*S is a diagonal matrix; + +// 1) +A=[1 2 -2 ;1 1 1;1 3 -1]; + +B=[1 0 0;0 1 0; 0 0 1]; + + [x,lam] = geigenvectors(A,B); + + inv(x)*A*x + + // 2) +A=[3 2 2;2 5 2;2 2 3]; + +B=[1 0 0;0 1 0; 0 0 1]; + + + [x,lam] = geigenvectors(A,B); + + inv(x)*A*x + + \ No newline at end of file diff --git a/50/CH3/EX3.4/ex_3_4.sce b/50/CH3/EX3.4/ex_3_4.sce new file mode 100755 index 000000000..9cbb29d89 --- /dev/null +++ b/50/CH3/EX3.4/ex_3_4.sce @@ -0,0 +1,30 @@ +// example no.3.4 +// solve the system of equations + +//(a) . by using cramer's rule, + +A=[1 2 -1;3 6 1;3 3 2] + +B1=[2 2 -1;1 6 1;3 3 2] + +B2=[1 2 -1;3 1 1;3 3 2] + +B3=[1 2 2;3 6 1; 3 3 3] + + +//we know; + +X1=det(B1)/det(A) +X2=det(B2)/det(A) +X3=det(B3)/det(A) + + +//(b). by detemining the inverse of the coefficient matrix + +A=[1 2 -1;3 6 1;3 3 2] + +b=[2 ;1 ;3] + +//we know; + +X=inv(A)*b \ No newline at end of file diff --git a/50/CH3/EX3.5/ex_3_5.sce b/50/CH3/EX3.5/ex_3_5.sce new file mode 100755 index 000000000..7f23bc40e --- /dev/null +++ b/50/CH3/EX3.5/ex_3_5.sce @@ -0,0 +1,9 @@ +//example-3.5 +//caption-solution by gauss elimibnation method + +A=[10 -1 2;1 10 -1;2 3 20] //matrices A and b from the above + //system of equations +b=[4;3;7] + +gausselim(A,b) //call gauss elimination function to solve the + //matrices A and b \ No newline at end of file diff --git a/50/CH3/EX3.6/ex_3_6.sce b/50/CH3/EX3.6/ex_3_6.sce new file mode 100755 index 000000000..47ca03a96 --- /dev/null +++ b/50/CH3/EX3.6/ex_3_6.sce @@ -0,0 +1,10 @@ +//example-3.6 +//caption-solution by gauss elimibnation method + +A=[1 1 1;3 3 4;2 1 3] //matrices A and b from the above + //system of equations + +b=[6;20;13] + +pivotgausselim(A,b) //call gauss elimination function to solve the + //matrices A and b \ No newline at end of file diff --git a/50/CH3/EX3.8/ex_3_8.sce b/50/CH3/EX3.8/ex_3_8.sce new file mode 100755 index 000000000..88c13b1ea --- /dev/null +++ b/50/CH3/EX3.8/ex_3_8.sce @@ -0,0 +1,8 @@ +//example no. 3.8 +// solving the matrix equation with partial pivoting in gauss elimination + +A=[2 1 1 -2;4 0 2 1;3 2 2 0;1 3 2 -1] + +b=[-10;8;7;-5] + +pivotgausselim(A,b) \ No newline at end of file diff --git a/50/CH3/EX3.9/ex_3_9.sce b/50/CH3/EX3.9/ex_3_9.sce new file mode 100755 index 000000000..34c7609a3 --- /dev/null +++ b/50/CH3/EX3.9/ex_3_9.sce @@ -0,0 +1,15 @@ +//example no.3.9 +//solving the system using inverse of the cofficient matrix + +A=[1 1 1;4 3 -1;3 5 3] + +I=[1 0 0;0 1 0;0 0 1] + +b=[1 ;6 ;4] + +M=jorden(A,I) + +IA=M(1:3,4:6) + +X=IA*b + diff --git a/50/CH4/EX4.15/ex_4_15.sce b/50/CH4/EX4.15/ex_4_15.sce new file mode 100755 index 000000000..be4304926 --- /dev/null +++ b/50/CH4/EX4.15/ex_4_15.sce @@ -0,0 +1,26 @@ +// example 4.15: +// obtain the interpolate using backward differences polinomial + + +xL=[.1 .2 .3 .4 .5 ]' + +f=[1.4 1.56 1.76 2 2.28]' +n=2; + + +// hence; +disp('P=1.4+(x-.1)*(.16/.1)+(x-.5)*(x-.4)*(.04/.02)') +disp('P=2x^2+x+1.28'); + +// 1) obtain the interpolate at x=0.25; +x=0.25; +[P]=NBDP(x,n,xL,f); +P +disp('f(.25)=1.655'); + + +// 2) obtain the interpolate at x=0.35; +x=0.35; +[P]=NBDP(x,n,xL,f); +P +disp('f(.35)=1.875'); \ No newline at end of file diff --git a/50/CH4/EX4.20/ex_4_20.sce b/50/CH4/EX4.20/ex_4_20.sce new file mode 100755 index 000000000..77cba39f6 --- /dev/null +++ b/50/CH4/EX4.20/ex_4_20.sce @@ -0,0 +1,14 @@ +// example 4.20; +// hermite interpolation: + +x=[-1 0 1]; + +f=[1 1 3]; + +fp=[-5 1 7]; + + P= hermiteinterpol(x,f,fp); + +// hence; +disp('f(-0.5)=3/8'); +disp('f(0.5)=11/8'); \ No newline at end of file diff --git a/50/CH4/EX4.21/ex_4_21.sce b/50/CH4/EX4.21/ex_4_21.sce new file mode 100755 index 000000000..288cd93e5 --- /dev/null +++ b/50/CH4/EX4.21/ex_4_21.sce @@ -0,0 +1,17 @@ +// example: 4.21; +// piecewise linear interpolating polinomials: + +x1=1;x2=2; x3=4;x4=8; +f1=3; f2=7; f3=21; f4=73; +// we need to apply legranges interpolation in sub-ranges [1,2];[2,4],[4,8]; + + x=poly(0,"x"); + +P1=legrangeinterpol (x1,x2,f1,f2); // in the range [1,2] +P1 + +P1=legrangeinterpol (x2,x3,f2,f3); // in the range [2,4] +P1 + +P1=legrangeinterpol (x3,x4,f3,f4); // in the range [4,8] +P1 \ No newline at end of file diff --git a/50/CH4/EX4.22/ex_4_22.sce b/50/CH4/EX4.22/ex_4_22.sce new file mode 100755 index 000000000..972fd5058 --- /dev/null +++ b/50/CH4/EX4.22/ex_4_22.sce @@ -0,0 +1,30 @@ +// example: 4.22; +// piecewise quadratic interpolating polinomials: + +X=[-3 -2 -1 1 3 6 7]; +F=[369 222 171 165 207 990 1779]; +// we need to apply legranges interpolation in sub-ranges [-3 ,-1];[-1,3],[3,7]; + + x=poly(0,"x"); + + // 1) in the range [-3,-1] + x=[-3 -2 -1]; + f=[369 222 171]; + n=2; +P2=lagrangefundamentalpoly(x,f,n); + + // 2) in the range [-1,3] + x=[-1 1 3]; + f=[171 165 207]; +n=2; +P2=lagrangefundamentalpoly(x,f,n) + + // 3) in the range [3,7] + x=[3 6 7]; + f=[207 990 1779]; +n=2; +P2=lagrangefundamentalpoly(x,f,n) + + + +// hence, we obtain the values of f(-2.5)=48; f(6.5)=1351.5; \ No newline at end of file diff --git a/50/CH4/EX4.23/ex_4_23.sce b/50/CH4/EX4.23/ex_4_23.sce new file mode 100755 index 000000000..e06601bef --- /dev/null +++ b/50/CH4/EX4.23/ex_4_23.sce @@ -0,0 +1,26 @@ +// example: 4.23; +// piecewise cubical interpolating polinomials: + +X=[-3 -2 -1 1 3 6 7]; +F=[369 222 171 165 207 990 1779]; +// we need to apply legranges interpolation in sub-ranges [-3 ,1];[1,7]; + + + x=poly(0,"x"); + + // 1) in the range [-3,1] + x=[-3 -2 -1 1]; + f=[369 222 171 165]; + n=3; +P2=lagrangefundamentalpoly(x,f,n); + + // 2) in the range [1,7] + x=[1 3 6 7]; + f=[165 207 990 1779]; +n=3; +P2=lagrangefundamentalpoly(x,f,n) + + + + // hence, + disp('f(6.5)=1339.25'); \ No newline at end of file diff --git a/50/CH4/EX4.3/ex_4_3.sce b/50/CH4/EX4.3/ex_4_3.sce new file mode 100755 index 000000000..fb99fe436 --- /dev/null +++ b/50/CH4/EX4.3/ex_4_3.sce @@ -0,0 +1,19 @@ + // example 4.3 + // find the linear interpolation polinomial + // using + + disp('f(2)=4'); + disp('f(2.5)=5.5'); + // 1)lagrange interpolation, + + P1=legrangeinterpol (2,2.5,4,5.5) + // 2)aitken's iterated interpolation, + + P1=aitkeninterpol (2,2.5,4,5.5) + + // 3) newton devided differance interpolation, + + P1=NDDinterpol (2,2.5,4,5.5) + + // hence approximate value of f(2.2)= 4.6; + \ No newline at end of file diff --git a/50/CH4/EX4.31/ex_4_31.sce b/50/CH4/EX4.31/ex_4_31.sce new file mode 100755 index 000000000..117c4dfb9 --- /dev/null +++ b/50/CH4/EX4.31/ex_4_31.sce @@ -0,0 +1,23 @@ +// example 4.31 +// obtain the linear polinomial approximation to the function f(x)=x^3 + +// let P(x)=a0*x+a1 + +// hence I(a0,a1)= integral (x^3-(a0*x+a1))^2 in the interval [0,1] + +printf('I=1/7-2*(a0/5+a1/4)+a0^2/3+a0*a1+a1^2') +printf('dI/da0 = -2/5+2/3*a0+a1=0') + +printf('dI/da1 = -1/2+a0+2*a1=0') + +// hence + +printf('[2/3 1;1 2]*[a0 ;a1]=[2/5; 1/2]') + +// solving for a0 and a1; + +a0=9/10; +a1=-1/5; +// hence considering the polinomial with intercept P1(x)=(9*x-2)/10; + +// considering the polinomial approximation through origin P2(x)=3*x/5; \ No newline at end of file diff --git a/50/CH4/EX4.32/ex_4_32.sce b/50/CH4/EX4.32/ex_4_32.sce new file mode 100755 index 000000000..81c43db58 --- /dev/null +++ b/50/CH4/EX4.32/ex_4_32.sce @@ -0,0 +1,45 @@ +// example 4.32 +// obtain the linear polinomial approximation to the function f(x)=x^1/2 + +// let P(x)=a0*x+a1 + + +// for n=1; +// hence I(c0,c1)= integral (x^1/2-(c1*x+c0))^2 in the interval [0,1] + + +printf('dI/dc0 = -2*(2/3-c0-c1/2)=0') + +printf('dI/dc1 =-2*(2/5-c0/2-c1/3) =0') + +// hence + +printf('[1 1/2;1/2 1/3]*[c0 ;c1]=[-4/3; -4/5]') + +// hence solving for c0 and c1; + + +// the first degree square approximation P(x)=4*(1+3*x)/15; + +// for n=2; + +// hence I(c0,c1,c2)= integral (x^1/2-(c2*x^2+c1*x+c0))^2 in the interval [0,1] + + +printf('dI/dc0 = (2/3-c0-c1/2-c2/2)=0') + +printf('dI/dc1 =(2/5-c0/2-c1/3-c2/4) =0') + +printf('dI/dc2 =(2/7-c0/3-c1/4-c2/5) =0') + + +// hence + +printf('[1 1/2 1/2;1/2 1/3 1/4;1/3 1/4 1/5]*[c0 ;c1;c2]=[-2/3; -2/5;-2/7]') + +// hence solving for c0,c1 and c2; + + +// the first degree square approximation P(x)=(6+48*x-20*x^2)/35; + + diff --git a/50/CH4/EX4.34/ex_4_34.sce b/50/CH4/EX4.34/ex_4_34.sce new file mode 100755 index 000000000..ed9752978 --- /dev/null +++ b/50/CH4/EX4.34/ex_4_34.sce @@ -0,0 +1,9 @@ +// example 4.34 + +// obtain least square straight line fit + +x=[.2 .4 .6 .8 1]; +f=[.447 .632 .775 .894 1]; + + +[P]=straightlineapprox(x,f) // call of the function to get the desired solution \ No newline at end of file diff --git a/50/CH4/EX4.35/ex_4_35.sce b/50/CH4/EX4.35/ex_4_35.sce new file mode 100755 index 000000000..d7dceb169 --- /dev/null +++ b/50/CH4/EX4.35/ex_4_35.sce @@ -0,0 +1,7 @@ +// example 4.35 + +// obtain least square approximation of second degree; +x=[-2 -1 0 1 2]; +f=[15 1 1 3 19]; + +[P]=quadraticapprox(x,f) // call of the function to get the desired solution \ No newline at end of file diff --git a/50/CH4/EX4.36/ex_4_36.sce b/50/CH4/EX4.36/ex_4_36.sce new file mode 100755 index 000000000..afc372683 --- /dev/null +++ b/50/CH4/EX4.36/ex_4_36.sce @@ -0,0 +1,31 @@ +// example 4.36 +// method of least squares to fit the data to the curve P(x)=c0*X+c1/squrt(X) + +x=[.2 .3 .5 1 2]; +f=[16 14 11 6 3]; + +// I(c0,c1)= summation of (f(x)-(c0*X+c1/sqrt(X)))^2 + +// hence on parcially derivating the summation, + +n=length(x);m=length(f); +if m<>n then + error('linreg - Vectors x and f are not of the same length.'); + abort; +end; + +s1=0; // s1= summation of x(i)*f(i) +s2=0; // s2= summation of f(i)/sqrt(x(i)) +s3=0; +for i=1:n + s1=s1+x(i)*f(i); + s2=s2+f(i)/sqrt(x(i)); + s3=s3+1/x(i); +end + +c0=det([s1 sum(sqrt(x));s2 s3])/det([sum(x^2) sum(sqrt(x));sum(sqrt(x)) s3]) + +c1=det([sum(x^2) s1;sum(sqrt(x)) s2])/det([sum(x^2) sum(sqrt(x));sum(sqrt(x)) s3]) +X=poly(0,"X"); +P=c0*X+c1/X^1/2 +// hence considering the polinomial P(x)=7.5961*X^1/2-1.1836*X \ No newline at end of file diff --git a/50/CH4/EX4.37/ex_4_37.sce b/50/CH4/EX4.37/ex_4_37.sce new file mode 100755 index 000000000..067c73c44 --- /dev/null +++ b/50/CH4/EX4.37/ex_4_37.sce @@ -0,0 +1,30 @@ +// example 4.37 +// method of least squares to fit the data to the curve P(x)=a*%e^(-3*t)+b*%e^(-2*t); + +t=[.1 .2 .3 .4]; +f=[.76 .58 .44 .35]; + +// I(c0,c1)= summation of (f(x)-a*%e^(-3*t)+b*%e^(-2*t)) + +// hence on parcially derivating the summation, + +n=length(t);m=length(f); +if m<>n then + error('linreg - Vectors t and f are not of the same length.'); + abort; +end; + +s1=0; // s1= summation of f(i)*%e^(-3*t(i)); +s2=0; // s2= summation of f(i)*%e^(-2*t(i)); + +for i=1:n + s1=s1+f(i)*%e^(-3*t(i)); + s2=s2+f(i)*%e^(-2*t(i)); + +end + +a=det([s1 sum(%e^(-5*t)); s2 sum(%e^(-4*t))])/det([sum(%e^(-6*t)) sum(%e^(-5*t)); sum(%e^(-5*t)) sum(%e^(-4*t))]) + +b=det([sum(%e^(-6*t)) s1; sum(%e^(-5*t)) s2])/det([sum(%e^(-6*t)) sum(%e^(-5*t)); sum(%e^(-5*t)) sum(%e^(-4*t))]) + +// hence considering the polinomial P(t)=.06853*%e^(-3*t)+0.3058*%e^(-2*t) \ No newline at end of file diff --git a/50/CH4/EX4.38/ex_4_38.sce b/50/CH4/EX4.38/ex_4_38.sce new file mode 100755 index 000000000..96d5857a9 --- /dev/null +++ b/50/CH4/EX4.38/ex_4_38.sce @@ -0,0 +1,26 @@ +// example 4.38 +// gram schmidt orthogonalisation + +W=1; + x=poly(0,"x"); + P0=1; + phi0=P0; + a10=integrate('W*x*phi0','x',0,1)/integrate('W*1*phi0','x',0,1) + P1=x-a10*phi0 + phi1=P1; + + a20=integrate('W*x^2*phi0','x',0,1)/integrate('W*1*phi0','x',0,1) + + a21=integrate('(x^2)*(x-1/2)','x',0,1)/integrate('(x-1/2)^2','x',0,1) + + P2=x^2-a20*x-a21*phi1 + +// since ,I= intgral [x^(1/2)-c0*P0-c1*P1-c2*P2]^2 inthe range [0,1] + +// hence partially derivating I + +c0=integrate('x^(1/2)','x',0,1)/integrate('1','x',0,1) +c1=integrate('(x^(1/2))*(x-(1/2))','x',0,1)/integrate('(x-(1/2))^2','x',0,1) +c1=integrate('(x^(1/2))*(x^2-4*x/3+1/2)','x',0,1)/integrate('(x^2-4*x/3+1/2)^2','x',0,1) + + diff --git a/50/CH4/EX4.39/ex_4_39.sce b/50/CH4/EX4.39/ex_4_39.sce new file mode 100755 index 000000000..3b231135d --- /dev/null +++ b/50/CH4/EX4.39/ex_4_39.sce @@ -0,0 +1,34 @@ +// example 4.39 +// gram schmidt orthogonalisation + +// 1) + W=1; + x=poly(0,"x"); + P0=1 + phi0=P0; + a10=0; + P1=x-a10*phi0 + phi1=P1; + + a20=integrate('x^2','x',-1,1)/integrate('W*1*phi0','x',-1,1); + + a21=integrate('(x^3)','x',-1,1)/integrate('(x)^2','x',-1,1); + + P2=x^2-a20*x-a21*phi1 + + +// 2) +disp(' W=1/(1-x^2)^(1/2)'); + x=poly(0,"x"); + P0=1 + phi0=P0; + a10=0; + P1=x-a10*phi0 + phi1=P1; + + a20=integrate('x^2/(1-x^2)^(1/2)','x',-1,1)/integrate('1/(1-x^2)^(1/2)','x',-1,1); + + a21=0; // since x^3 is an odd function; + + P2=x^2-a20*x-a21*phi1 + diff --git a/50/CH4/EX4.4/ex_4_4.sce b/50/CH4/EX4.4/ex_4_4.sce new file mode 100755 index 000000000..c31365e30 --- /dev/null +++ b/50/CH4/EX4.4/ex_4_4.sce @@ -0,0 +1,12 @@ + // example 4.4 + // find the linear interpolation polinomial + // using lagrange interpolation, + + disp('sin(.1)=.09983; sin(.2)=.19867'); + + + P1=legrangeinterpol (.1,.2,.09983,.19867) + +// hence; +disp('P(.15)=.00099+.9884*.15') +disp('P(0.15)=0.14925'); \ No newline at end of file diff --git a/50/CH4/EX4.41/ex_4_41.sce b/50/CH4/EX4.41/ex_4_41.sce new file mode 100755 index 000000000..1ce2864ac --- /dev/null +++ b/50/CH4/EX4.41/ex_4_41.sce @@ -0,0 +1,20 @@ +// example 4.41 +// using chebyshev polinomials obtain least squares approximation of second degree; + +// the chebeshev polinomials; +x=poly(0,"x"); +T0=1; +T1=x; +T2=2*x^2-1; + + +// I=integrate ('1/(1-x^2)^(1/2)*(x^4-c0*T0-c1*T1-c2*T2)^2','x',-1,1) + +// since; +c0=integrate('(1/3.14)*(x^4)/(1-x^2)^(1/2)','x',-1,1) + +c1=integrate('(2/3.14)*(x^5)/(1-x^2)^(1/2)','x',-1,1) + +c2=integrate('(2/3.14)*(x^4)*(2*x^2-1)/(1-x^2)^(1/2)','x',-1,1) + + f=(3/8)*T0+(1/2)*T2; \ No newline at end of file diff --git a/50/CH4/EX4.6/ex_4_6.sce b/50/CH4/EX4.6/ex_4_6.sce new file mode 100755 index 000000000..869b1f772 --- /dev/null +++ b/50/CH4/EX4.6/ex_4_6.sce @@ -0,0 +1,18 @@ +// example :4.6 +// caption : obtain the legrange linear interpolating polinomial + + +// 1) obtain the legrange linear interpolating polinomial in the interval [1,3] and obtain approximate value of f(1.5),f(2.5); +x0=1;x1=2;x2=3;f0=.8415;f1=.9093;f2=.1411; + +P13=legrangeinterpol (x0,x2,f0,f2) //in the range [1 ,3] + + +P12=legrangeinterpol (x0,x1,f0,f1) // in the range [1,2] + +P23=legrangeinterpol (x1,x2,f1,f2) // inthe range [2,3] + +// from P23 we find that ; where as exact value is sin(2.5)=0.5985; +disp('P(1.5)=0.8754'); +disp('exact value of sin(1.5)=.9975') +disp('P(2.5)=0.5252'); \ No newline at end of file diff --git a/50/CH4/EX4.7/ex_4_7.sce b/50/CH4/EX4.7/ex_4_7.sce new file mode 100755 index 000000000..cbc217059 --- /dev/null +++ b/50/CH4/EX4.7/ex_4_7.sce @@ -0,0 +1,15 @@ + // example 4.7 + // polinomial of degree 2; + + // f(0)=1;f(1)=3;f(3)=55; + + // using legrange fundamental polinomial rule, + + x=[0 1 3]; // arrainging the inputs of the function as elements of a row, + f=[1 3 55]; // arrainging the outputs of the function as elements of a row, + n=2; // degree of the polinomial; + + + P2=lagrangefundamentalpoly(x,f,n) + + \ No newline at end of file diff --git a/50/CH4/EX4.8/ex_4_8.sce b/50/CH4/EX4.8/ex_4_8.sce new file mode 100755 index 000000000..6bfdd8e78 --- /dev/null +++ b/50/CH4/EX4.8/ex_4_8.sce @@ -0,0 +1,15 @@ + + // example 4.8 + // caption: solution by quadratic interpolation; + + // x-degrees:[10 20 30] + // hence x in radians is + x=[3.14/18 3.14/9 3.14/6]; + f=[1.1585 1.2817 1.3660]; + n=2; + + + P2=lagrangefundamentalpoly(x,f,n) + + // hence from P2 ,the exact value of f(3.14/12) is 1.2246; + // where as exact value is 1.2247; \ No newline at end of file diff --git a/50/CH4/EX4.9/ex_4_9.sce b/50/CH4/EX4.9/ex_4_9.sce new file mode 100755 index 000000000..726da7ed5 --- /dev/null +++ b/50/CH4/EX4.9/ex_4_9.sce @@ -0,0 +1,16 @@ +// example 4.9 +// caption :obtain the polinomial of degree 2 + +x=[0 1 3]; +f=[1 3 55]; +n=2; + +// 1) iterated interpolation; + + +[L012,L02,L01]=iteratedinterpol (x,f,n) + +// 2) newton divided diffrences interpolation; + + + P2=NDDinterpol2 (x,f) \ No newline at end of file diff --git a/50/CH5/EX5.1/ex_5_1.sce b/50/CH5/EX5.1/ex_5_1.sce new file mode 100755 index 000000000..4ffe53138 --- /dev/null +++ b/50/CH5/EX5.1/ex_5_1.sce @@ -0,0 +1,12 @@ +// example: 5.1 +// linear and quadratic interpolation: + +// f(x)=ln x; + +xL=[2 2.2 2.6]; +f=[.69315 .78846 .95551]; + +// 1) fp(2) with linear interpolation; + +fp=linearinterpol(xL,f); +disp(fp); diff --git a/50/CH5/EX5.10/ex_5_10.sce b/50/CH5/EX5.10/ex_5_10.sce new file mode 100755 index 000000000..919949c50 --- /dev/null +++ b/50/CH5/EX5.10/ex_5_10.sce @@ -0,0 +1,15 @@ +// example 5.10; +// find the jacobian matrix; + + +// given two functions in x,y; +// and the point at which the jacobian has to be found out; + +deff('[w]=f1(x,y)','w=x^2+y^2-x'); + +deff('[q]=f2(x,y)','q=x^2-y^2-y'); + +h=1;k=1; + +J= jacobianmat (f1,f2,h,k); +disp(J); \ No newline at end of file diff --git a/50/CH5/EX5.11/ex_5_11.sce b/50/CH5/EX5.11/ex_5_11.sce new file mode 100755 index 000000000..3f861a3ea --- /dev/null +++ b/50/CH5/EX5.11/ex_5_11.sce @@ -0,0 +1,18 @@ +// example : 5.11 +// solve the definite integral by 1) trapezoidal rule, 2)simpsons rule +// exact value of the integral is ln 2= 0.693147, + +deff('[y]=F(x)','y=1/(1+x)') + +//1)trapezoidal rule, + +a=0; +b=1; +I =trapezoidal(0,1,F) +disp(error =.75-.693147) + +// simpson's rule + +I=simpson(a,b,F) + +disp(error =.694444-.693147) \ No newline at end of file diff --git a/50/CH5/EX5.12/ex_5_12.sce b/50/CH5/EX5.12/ex_5_12.sce new file mode 100755 index 000000000..2e08e16f7 --- /dev/null +++ b/50/CH5/EX5.12/ex_5_12.sce @@ -0,0 +1,17 @@ +// example 5.12 +// caption: solve the integral by 1)mid-point rule,2)two-point open type rule + + +// let integration of f(x)=sin(x)/(x) in the range [0,1] is equal to I1 and I2 +// 1)mid -point rule; +a=0;b=1; +h=(b-a)/2; + +x=0:h:1; +deff('[y]=f(x)','y=sin(x)/x') +I1=2*h*f(x(1)+h) + + +//2) two-point open type rule +h=(b-a)/3; +I2=(3/2)*h*(f(x(1)+h)+f(x(1)+2*h)) \ No newline at end of file diff --git a/50/CH5/EX5.13/ex_5_13.sce b/50/CH5/EX5.13/ex_5_13.sce new file mode 100755 index 000000000..d0ada71c8 --- /dev/null +++ b/50/CH5/EX5.13/ex_5_13.sce @@ -0,0 +1,12 @@ +// example 5.13 +// caption: simpson 3-8 rule + + +// let integration of f(x)=1/(1+x) in the range [0,1] by simpson 3-8 rule is equal to I + +x=0:1/3:1; +deff('[y]=f(x)','y=1/(1+x)') + +[I] = simpson38(x,f) + + \ No newline at end of file diff --git a/50/CH5/EX5.15/ex_5_15.sce b/50/CH5/EX5.15/ex_5_15.sce new file mode 100755 index 000000000..b56793e06 --- /dev/null +++ b/50/CH5/EX5.15/ex_5_15.sce @@ -0,0 +1,26 @@ +// example :5.15 +// find the quadrature formula of +// integral of f(x)*(1/sqrt(x(1+x))) in the range [0,1]= a1*f(0)+a2*f(1/2)+a3*f(1)=I +// hence find integral 1/sqrt(x-x^3) in the range [0,1] + +// making the method exact for polinomials of degree upto 2, +// I=I1=a1+a2+a3 +// I=I2=(1/2)*a2+a3 +// I=I3=(1/4)*a2+a3 + +// A=[a1 a2 a3]' + +I1=integrate('1/sqrt(x*(1-x))','x',0,1) +I2=integrate('x/sqrt(x*(1-x))','x',0,1) +I3=integrate('x^2/sqrt(x*(1-x))','x',0,1) + +//hence +// [1 1 1;0 1/2 1 ;0 1/4 1]*A=[I1 I2 I3]' + +A=inv([1 1 1;0 1/2 1 ;0 1/4 1])*[I1 I2 I3]' +// I=(3.14/4)*(f(0)+2*f(1/2)+f(1)); + +// hence, for solving the integral 1/sqrt(x-x^3) in the range [0,1]=I + +deff('[y]=f(x)','y=1/sqrt(1+x)'); +I=(3.14/4)*[1+2*sqrt(2/3)+sqrt(2)/2] \ No newline at end of file diff --git a/50/CH5/EX5.16/ex_5_16.sce b/50/CH5/EX5.16/ex_5_16.sce new file mode 100755 index 000000000..f3b260f89 --- /dev/null +++ b/50/CH5/EX5.16/ex_5_16.sce @@ -0,0 +1,16 @@ +// example 5.16 +// caption: gauss-legendre three point method +// I= integral 1/(1+x) in the range [0,1]; +// first we need ti transform the interval [0,1 ] to [-1,1], since gauss-legendre three point method is applicable in the range[-1,1], + +// let t=ax+b; +// solving for a,b from the two ranges, we get a=2; b=-1; t=2x-1; + +// hence I=integral 1/(1+x) in the range [0,1]= integral 1/(t+3) in the range [-1,1]; + + +deff('[y]=f(t)','y=1/(t+3)'); +// since , from gauss legendre three point rule(n=2); +I=(1/9)*(5*f(-sqrt(3/5))+8*f(0)+5*f(sqrt(3/5))) + +// we know , exact solution is ln 2=0.693147; \ No newline at end of file diff --git a/50/CH5/EX5.17/ex_5_17.sce b/50/CH5/EX5.17/ex_5_17.sce new file mode 100755 index 000000000..f0e6d0b6f --- /dev/null +++ b/50/CH5/EX5.17/ex_5_17.sce @@ -0,0 +1,24 @@ +// example 5.17 +// caption: gauss-legendre method +// I= integral 2*x/(1+x^4) in the range [1,2]; +// first we need ti transform the interval [1,2 ] to [-1,1], since gauss-legendre three point method is applicable in the range[-1,1], + +// let t=ax+b; +// solving for a,b from the two ranges, we get a=1/2; b=3/2; x=(t+3)/2; + +// hence I=integral 2*x/(1+x^4) in the range [0,1]= integral 8*(t+3)/16+(t+3)^4 in the range [-1,1]; + + +deff('[y]=f(t)','y=8*(t+3)/(16+(t+3)^4) '); + +// 1) since , from gauss legendre one point rule; +I1=2*f(0) + +// 2) since , from gauss legendre two point rule; +I2=f(-1/sqrt(3))+f(1/sqrt(3)) + +// 3) since , from gauss legendre three point rule; +I=(1/9)*(5*f(-sqrt(3/5))+8*f(0)+5*f(sqrt(3/5))) + + +// we know , exact solution is 0.5404; \ No newline at end of file diff --git a/50/CH5/EX5.18/ex_5_18.sce b/50/CH5/EX5.18/ex_5_18.sce new file mode 100755 index 000000000..08ce9df7a --- /dev/null +++ b/50/CH5/EX5.18/ex_5_18.sce @@ -0,0 +1,23 @@ +// example 5.18 +// caption: gauss-chebyshev method + +// we write the integral as I=integral f(x)/sqrt(1-x^2) in the range [-1,1]; +// where f(x)=(1-x^2)^2*cos(x) + +deff('[y]=f(x)','y=(1-x^2)^2*cos(x)'); + +// 1) since , from gauss chebyshev one point rule; +I1=(3.14)*f(0) + +// 2) since , from gauss chebyshev two point rule; +I2=(3.14/2)*f(-1/sqrt(2))+f(1/sqrt(2)) + +// 3) since , from gauss chebyshev three point rule; +I=(3.14/3)*(f(-sqrt(3)/2)+f(0)+f(sqrt(3)/2)) + + +// and 4) since , from gauss legendre three point rule; +I=(1/9)*(5*f(-sqrt(3/5))+8*f(0)+5*f(sqrt(3/5))) + + + diff --git a/50/CH5/EX5.2/ex_5_2.sce b/50/CH5/EX5.2/ex_5_2.sce new file mode 100755 index 000000000..cd5e0032e --- /dev/null +++ b/50/CH5/EX5.2/ex_5_2.sce @@ -0,0 +1,9 @@ +// example 5.2 +// evaluate fp(.8) and fpp(.8) with quadratic interpolation; + +xL=[.4 .6 .8]; +f=[.0256 .1296 .4096]; +h=.2; + +fp=(1/2*h)*(f(1)-4*f(2)+3*f(3)) +fpp=(1/h^2)*(f(1)-2*f(2)+f(3)) // from equation 5.22c and 5.24c in the book; diff --git a/50/CH5/EX5.20/ex_5_20.sce b/50/CH5/EX5.20/ex_5_20.sce new file mode 100755 index 000000000..3e0f84277 --- /dev/null +++ b/50/CH5/EX5.20/ex_5_20.sce @@ -0,0 +1,16 @@ +// example 5.20 +// caption: gauss-leguerre method +// I= integral e^-x/(1+x^2) in the range [0,~]; + +// observing the integral we can inffer that f(x)=1/(1+x^2) + +deff('[y]=f(x)','y=1/(1+x^2) '); + + +// 1) since , from gauss leguerre two point rule; +I2=(1/4)*[(2+sqrt(2))*f(2-sqrt(2))+(2-sqrt(2))*f(2+sqrt(2))] + +// 3) since , from gauss leguerre three point rule; +I=(0.71109*f(0.41577)+0.27852*f(2.29428)+0.01039*f(6.28995)) + + diff --git a/50/CH5/EX5.21/ex_5_21.sce b/50/CH5/EX5.21/ex_5_21.sce new file mode 100755 index 000000000..64303a5b5 --- /dev/null +++ b/50/CH5/EX5.21/ex_5_21.sce @@ -0,0 +1,16 @@ +// example 5.21 +// caption: gauss-leguerre method +// I= integral e^-x*(3*x^3-5*x+1) in the range [0,~]; + +// observing the integral we can inffer that f(x)=(3*x^3-5*x+1) + +deff('[y]=f(x)','y=(3*x^3-5*x+1) '); + + +// 1) since , from gauss leguerre two point rule; +I2=(1/4)*[(2+sqrt(2))*f(2-sqrt(2))+(2-sqrt(2))*f(2+sqrt(2))] + +// 3) since , from gauss leguerre three point rule; +I3=(0.71109*f(0.41577)+0.27852*f(2.29428)+0.01039*f(6.28995)) + + diff --git a/50/CH5/EX5.22/ex_5_22.sce b/50/CH5/EX5.22/ex_5_22.sce new file mode 100755 index 000000000..7afdb7036 --- /dev/null +++ b/50/CH5/EX5.22/ex_5_22.sce @@ -0,0 +1,20 @@ +// example 5.22 +// caption: gauss-leguerre method +// I= integral 1/(x^2+2*x+2) in the range [0,~]; + +// since in the gauss-leguerre method the integral would be of the form e^x*f(x); + +// observing the integral we can inffer that f(x)=%e^x/(x^2+2*x+2) +deff('[y]=f(x)','y=%e^x/(x^2+2*x+2) '); + + +// 1) since , from gauss leguerre two point rule; +I2=(1/4)*[(2+sqrt(2))*f(2-sqrt(2))+(2-sqrt(2))*f(2+sqrt(2))] + +// 3) since , from gauss leguerre three point rule; +I=(0.71109*f(0.41577)+0.27852*f(2.29428)+0.01039*f(6.28995)) + + +// the exact solution is given by, + +I=integrate('1/((x+1)^2+1)','x',0,1000) // 1000 ~infinite; diff --git a/50/CH5/EX5.26/ex_5_26.sce b/50/CH5/EX5.26/ex_5_26.sce new file mode 100755 index 000000000..402f55cb0 --- /dev/null +++ b/50/CH5/EX5.26/ex_5_26.sce @@ -0,0 +1,42 @@ +// Example 5.26 +// caption: 1) composite trapizoidal rule, 2) composite simpsons rule with 2,4 ,8 equal sub-intervals, + +// I=integral 1/(1+x) in the range [0,1] + +deff('[y]=f(x)','y=1/(1+x)') + +// when N=2; +// 1)composite trapizoidal rule +h=1/2; +x=0:h:1; + + IT=comptrapezoidal(x,h,f) + + // 2)composite simpsons rule + + [I] = simpson13(x,h,f) + + + // when N=4 + // 1)composite trapizoidal rule +h=1/4; +x=0:h:1; + + IT=comptrapezoidal(x,h,f) + + // 2)composite simpsons rule + + [I] = simpson13(x,h,f) + + + + // when N=8 + // 1)composite trapizoidal rule +h=1/8; +x=0:h:1; + + IT=comptrapezoidal(x,h,f) + + // 2)composite simpsons rule + + [I] = simpson13(x,h,f) \ No newline at end of file diff --git a/50/CH5/EX5.27/ex_5_27.sce b/50/CH5/EX5.27/ex_5_27.sce new file mode 100755 index 000000000..cd5800369 --- /dev/null +++ b/50/CH5/EX5.27/ex_5_27.sce @@ -0,0 +1,23 @@ +// example 5.27 +// caption: gauss-legendre three point method +// I= integral 1/(1+x) in the range [0,1]; + +// we are asked to subdivide the range into two, +// first we need to sub-divide the interval [0,1 ] to [0,1/2] and [1/2,1] and then transform both to [-1,1], since gauss-legendre three point method is applicable in the range[-1,1], + +// t=4x-1 and y=4x-3; + +// hence I=integral 1/(1+x) in the range [0,1]= integral 1/(t+5) in the range [-1,1]+ integral 1/(t+7) in the range [-1,1] + + +deff('[y1]=f1(t)','y1=1/(t+5)'); +// since , from gauss legendre three point rule(n=2); +I1=(1/9)*(5*f1(-sqrt(3/5))+8*f1(0)+5*f1(sqrt(3/5))) + +deff('[y2]=f2(t)','y2=1/(t+7)'); +// since , from gauss legendre three point rule(n=2); +I2=(1/9)*(5*f2(-sqrt(3/5))+8*f2(0)+5*f2(sqrt(3/5))) + +I=I1+I2 + +// we know , exact solution is .693147; \ No newline at end of file diff --git a/50/CH5/EX5.29/ex_5_29.sce b/50/CH5/EX5.29/ex_5_29.sce new file mode 100755 index 000000000..7e88c2c8c --- /dev/null +++ b/50/CH5/EX5.29/ex_5_29.sce @@ -0,0 +1,12 @@ +// example 5.29 +// evaluate the given double integral using the simpsons rule; + +// I= double integral f(x)=1/(x+y) in the range x=[1,2],y=[1,1.5]; + +h=.5; +k=.25; +deff('[w]=f(x,y)','w=1/(x+y)') + +I=(.125/9)*[{f(1,1)+f(2,1)+f(1,1.5)+f(2,1.5)}+4*{f(1.5,1)+f(1,1.25)+f(1.5,1.5)+f(2,1.25)}+16*f(1.5,1.25)]; +disp(I); + diff --git a/50/CH5/EX5.30/ex_5_30.sce b/50/CH5/EX5.30/ex_5_30.sce new file mode 100755 index 000000000..17f88a29d --- /dev/null +++ b/50/CH5/EX5.30/ex_5_30.sce @@ -0,0 +1,18 @@ +// example 5.30 +// evaluate the given double integral using the simpsons rule; + +// I= double integral f(x)=1/(x+y) in the range x=[1,2],y=[1,2]; +// 1) +h=.5; +k=.5; +deff('[w]=f(x,y)','w=1/(x+y)') + +I=(1/16)*[{f(1,1)+f(2,1)+f(1,2)+f(2,2)}+2*{f(1.5,1)+f(1,1.5)+f(2,1.5)+f(1.5,2)}+4*f(1.5,1.5)] + +// 2) +h=.25; +k=.25; +deff('[w]=f(x,y)','w=1/(x+y)') + +I=(1/64)*[{f(1,1)+f(2,1)+f(1,2)+f(2,2)}+2*{f(5/4,1)+f(3/2,1)+f(7/4,1)+f(1,5/4)+f(1,3/2)+f(1,7/4)+f(2,5/4)+f(2,3/4)+f(2,7/4)+f(5/4,2)+f(3/2,2)+f(7/4,2)}+4*{f(5/4,5/4)+f(5/4,3/2)+f(5/4,7/4)+f(3/2,5/4)+f(3/2,3/2)+f(3/2,7/4)+f(7/4,5/4)+f(7/4,3/2)+f(7/4,7/4)}] + diff --git a/50/CH6/EX6.12/ex_6_12.sce b/50/CH6/EX6.12/ex_6_12.sce new file mode 100755 index 000000000..045ec5415 --- /dev/null +++ b/50/CH6/EX6.12/ex_6_12.sce @@ -0,0 +1,9 @@ +// example 6.12, +// caption: solve the IVP by backward euler method, +// with h=0.2, +// u'=f(t,u) +deff('[z]=f(t,u)','z=-2*t*u^2'); + + +[u] = backeuler(1,0,0.4,.2,f) // h=0.2; + diff --git a/50/CH6/EX6.13/ex_6_13.sce b/50/CH6/EX6.13/ex_6_13.sce new file mode 100755 index 000000000..7fc402959 --- /dev/null +++ b/50/CH6/EX6.13/ex_6_13.sce @@ -0,0 +1,11 @@ +// example 6.13, +// caption: solve the IVP by euler midpoint method, +// with h=0.2, +// u'=f(t,u) +deff('[z]=f(t,u)','z=-2*t*u^2'); +deff('[w]=fp(t,u)','w=-2*u^2-4*u*t'); + + + +[u] = eulermidpoint(1,0,1,.2,f,fp) // h=0.2; + diff --git a/50/CH6/EX6.15/ex_6_15.sce b/50/CH6/EX6.15/ex_6_15.sce new file mode 100755 index 000000000..30daf63a0 --- /dev/null +++ b/50/CH6/EX6.15/ex_6_15.sce @@ -0,0 +1,20 @@ +// example 6.15 +//caption: solving ODE by tailor series method +// u'=t^2+ u^2, u0=0; + +t=0;U=0; //at t=0, the value of u is 0 +U1=0; // u1 is the 1st derivatove of the funtion u +U2=2*t+2*U*U1 // U2 ---- 2nd derivative +U3=2+2*(U*U2+U1^2) +U4=2*(U*U3+3*U1*U2) +U5=2*(U*U4+4*U1*U3+3*U2^2) +U6=2*(U*U5+5*U1*U4+10*U2*U3) +U7=2*(U*U6+6*U1*U5+15*U2*U4+10*U3^2) +U8=0; +U9=0; +U10=0; +U11=2*(U*U10+10*U1*U9+45*U2*U8+120*U3*U7+210*U4*U6+126*U5^2) + // U11 is the 11th derivative of u + + +taylor(1) \ No newline at end of file diff --git a/50/CH6/EX6.17/ex_6_17.sce b/50/CH6/EX6.17/ex_6_17.sce new file mode 100755 index 000000000..9de434d2a --- /dev/null +++ b/50/CH6/EX6.17/ex_6_17.sce @@ -0,0 +1,17 @@ +// example 6.17, +// caption: solve by 1)modified euler cauchy, 2)heun method + +// h=0.2 +// 1) modified euler cauchy method, + +// u'=f(t,u) +// u'=-2tu^2 +deff('[z]=f(t,u)','z=-2*t*u^2'); + +modifiedeuler(1,0,.4,.2,f) // calling the function, + +// 2) heun method, +deff('[z]=f(t,u)','z=-2*t*u^2'); + + +heun(1,0,.4,.2,f) // calling the function, diff --git a/50/CH6/EX6.18/ex_6_18.sce b/50/CH6/EX6.18/ex_6_18.sce new file mode 100755 index 000000000..5a8165e09 --- /dev/null +++ b/50/CH6/EX6.18/ex_6_18.sce @@ -0,0 +1,8 @@ +// example 6.18, +// caption: use of 4th order runge kutta method, + +// u'=f(t,u) +// u'=-2tu^2 +deff('[z]=f(t,u)','z=-2*t*u^2'); + +RK4(1,0,.4,.2,f) // calling the function, \ No newline at end of file diff --git a/50/CH6/EX6.20/ex_6_20.sce b/50/CH6/EX6.20/ex_6_20.sce new file mode 100755 index 000000000..5430115f1 --- /dev/null +++ b/50/CH6/EX6.20/ex_6_20.sce @@ -0,0 +1,18 @@ +// example no. 6.20, +// caption: solve the system of equations + +//1) eulercauchy method solving simultanious ODE + +deff('[z]=f1(t,u,v)','z=-3*u+2*v'); +deff('[w]=f2(t,u,v)','w=3*u-4*v'); + + +[u,v,t] = simeulercauchy(0,.5,0,.4,.2,f1,f2) + +// 2) RK4 method solving simultanious ODE + + + +[u,v,t]=simRK4(0,.5,0,.4,.2,f1,f2) + + diff --git a/50/CH6/EX6.21/ex_6_21.sce b/50/CH6/EX6.21/ex_6_21.sce new file mode 100755 index 000000000..32b9e0d95 --- /dev/null +++ b/50/CH6/EX6.21/ex_6_21.sce @@ -0,0 +1,27 @@ +// example no. 6.21, +// caption: solving the IVP by implicit RK2 method + +// u'=f(t,u) +// u'=-2tu^2 +//u(0)=1,h=0.2; +t0=0;h=0.2;tn=.4;u0=1; +deff('[z]=f(t,u)','z=-2*t*u^2'); +umaxAllowed = 1e+100; + +t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0; + +for j = 1:n-1 + // k1=h*f(t(j)+h/2,u(j)+k1/2); + // conidering the IVP we can infer that the above expression in non linear in k1, +// hence we use newton rapson method to solve for k1; +deff('[w]=F(u2)','w=k1+h*(2*t(j)+h)*(u(j)+k1/2)^2') // u2=u(2) +deff('[x]=Fp(u2)','x=1+h*(2*t(j)+h)*(u(j)+k1/2)') + +k1=h*f(t(j),u(j)); + +newton(k1,F,Fp); + u(j+1) = u(j) +k1 + disp(u(j+1)) + +end; + diff --git a/50/CH6/EX6.25/ex_6_25.sce b/50/CH6/EX6.25/ex_6_25.sce new file mode 100755 index 000000000..d9a405146 --- /dev/null +++ b/50/CH6/EX6.25/ex_6_25.sce @@ -0,0 +1,8 @@ +// example 6.25 +// caption: solving the IVP by adams-bashforth 3rd order method. +// u'=f(t,u) +// u'=-2tu^2 +//u(0)=1,h=0.2; +deff('[z]=f(t,u)','z=-2*t*u^2'); + +adamsbashforth3(1,0,1,.2,f) // calling the function, diff --git a/50/CH6/EX6.27/ex_6_27.sce b/50/CH6/EX6.27/ex_6_27.sce new file mode 100755 index 000000000..41081a0a9 --- /dev/null +++ b/50/CH6/EX6.27/ex_6_27.sce @@ -0,0 +1,32 @@ +// example 6.27 +// solving IVP by 3rd order adams moulton +// u'=t^2+u^2, u(1)=2, +// h=0.1, [1,1.2] +deff('[z]=f(t,u)','z=t^2+u^2'); +t0=1;u0=2;h=0.1;tn=1.2; +// third order adams moulton method, +// u(j+2)=u(j+1)+(h/12)*(5*f(t(j+2),u(j+2))+8*f(t(j+1),u(j+1))-f(t(j),u(j)));- is the expression for adamsbas-moulton3 + + +// on observing the IVP we can inffer that this would be a non linear equation, +// u(j+2)=u(j+1)+(h/12)*(5*((t(j+2))^2+(u(j+2))^2)+8*((t(j+1))^2+(u(j+1))^2)-((t(j))^2+(u(j))^2)) + +t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0; +for j = 1:n-2 + if j==1 then + k1=h*f(t(j),u(j)); + k2=h*f(t(j)+h,u(j)+k1); + u(j+1) = u(j) + (k2+k1)/2; + disp(u(j+1)) + end; +end; + +// hence the third order adams moulton expression turns to be, +// u(2)= 0.041667*(u(2))^2+3.194629 +// let us use newton raphsom method to solve this, +deff('[w]=F(u2)','w=-u2+ 0.041667*(u2)^2+3.194629') // u2=u(2) +deff('[x]=Fp(u2)','x=-1+ 0.041667*2*u2') + +// let us assume the initial guess of u(2)=u(1); + +newton(2.633333,F,Fp) \ No newline at end of file diff --git a/50/CH6/EX6.3/ex_6_3.sce b/50/CH6/EX6.3/ex_6_3.sce new file mode 100755 index 000000000..ea77769e4 --- /dev/null +++ b/50/CH6/EX6.3/ex_6_3.sce @@ -0,0 +1,17 @@ +// example 6.3 +// solution to the given IVP + +disp('du/dt= A*u'); +// u=[u1 u2]'; +A=[-3 4 ;-2 3]; // given +B=[1 0;0 1]; // identity matrix; + + + + +[x,lam] = geigenvectors(A,B); + +// hence; +disp('u=c1*%e^t*x(:,1)+c2*%e^-t*x(:,2)'); +disp('u1=c1*%e^t+c2*%e^-t*2') +disp('u2=c1*%e^t+c2*%e^-t') \ No newline at end of file diff --git a/50/CH6/EX6.32/ex_6_32.sce b/50/CH6/EX6.32/ex_6_32.sce new file mode 100755 index 000000000..eca7419c7 --- /dev/null +++ b/50/CH6/EX6.32/ex_6_32.sce @@ -0,0 +1,25 @@ +// example 6.32 +// caption: solving the IVP by numerov method +// u''=(1+t^2)*u +// u(0)=1, u'(0)=0 ,[0,1] +// h=0.2, + +// expression for numerov method is +//u(j+1)-2*u(j)+u(j-1)=(h^2/12)*(u''(j+1)+10*u''(j)+u''(j-1)); + +// observing the IVP we can reduce the numerov method to +//u(2)=2*u(1)-u(0)+(.2^2/12)*(1.16*u(2)+10.4*u(1)+1); for j=1 +// u(3)=2*u(2)-u(1)+(.2^2/12)*(1.36*u(3)+11.6*u(2)+1.04*u(1)); for j=2 +// u(4)=2*u(3)-u(2)+(.2^2/12)*(1.64*u(4)+13.6*u(3)+1.16*u(2)); for j=3 +// u(5)=2*u(4)-u(3)+(.2^2/12)*(2*u(5)+16.4*u(4)+1.36*u(3)); for j=4 + +// from taylor series expansion we observe that +u1=1.0202;u0=1; +//u2-(.2^2/12)*(1.16*u2)=2*u1-u0+(.2^2/12)*(10.4*u1+1); +u2=(1/.9961333)*2*u1-u0+(.2^2/12)*(10.4*u1+1) + +u3=(1/.995467)*(2.038667*u2-.996533*u1) + +u4=(1/.994533)*(2.045333*u3-.996133*u2) + +u5=(1/.993333)*(2.054667*u4-.995467*u3) \ No newline at end of file diff --git a/50/CH6/EX6.4/ex_6_4.sce b/50/CH6/EX6.4/ex_6_4.sce new file mode 100755 index 000000000..f24a42d57 --- /dev/null +++ b/50/CH6/EX6.4/ex_6_4.sce @@ -0,0 +1,16 @@ +// example 6.4 +// solution to the given IVP + +disp('du/dt= A*u'); +// u=[u1 u2]'; +A=[-2 1;1 -20]; // given +B=[1 0;0 1]; // identity matrix; + + + + + +[x,lam] = geigenvectors(A,B); + +// hence; +disp('u=c1*%e^(lam(1)*t)*x(:,1)+c2*%e^-(lam(2)*t)*x(:,2)'); diff --git a/50/CH6/EX6.9/ex_6_9.sce b/50/CH6/EX6.9/ex_6_9.sce new file mode 100755 index 000000000..b3b8ddf45 --- /dev/null +++ b/50/CH6/EX6.9/ex_6_9.sce @@ -0,0 +1,14 @@ +// example 6.9 +// solve the IVP by euler method, +// with h=0.2, 0.1, 0.05; +// u'=f(t,u) +// u'=-2tu^2 +deff('[z]=f(t,u)','z=-2*t*u^2'); + + +[u,t] = Euler1(1,0,1,.2,f) // h=0.2; + +[u,t] = Euler1(1,0,1,0.1,f) // h=0.1; + + +[u,t] = Euler1(1,0,1,0.05,f) // h=0.05; \ No newline at end of file diff --git a/50/CH7/EX7.1/ex_7_1.sce b/50/CH7/EX7.1/ex_7_1.sce new file mode 100755 index 000000000..d857f64d6 --- /dev/null +++ b/50/CH7/EX7.1/ex_7_1.sce @@ -0,0 +1,27 @@ +// example 7.1 +// solve by shooting method; + +// u''=u+1; +// u(0)=0; u(1)=%e-1; + +// let -> U1(x)=du/dx; +// U2(x)=d2u/dx2; + +// U(x)=[U1(x);U2(x)] + +// hence ; +// dU/dx=f(x,U); + + + +deff('[w]=f(x,U)','w=[U(2); U(1)+1]') + +h=0.25; +x=[0:h:1]; +ub=[0,%e-1]; +up=[0:1:10]; + + +[U] = shooting(ub,up,x,f); + +// the solution obtained would show the values of u and their derivatives at various x taken in regular intervals of h; \ No newline at end of file diff --git a/50/CH7/EX7.11/ex_7_11.sce b/50/CH7/EX7.11/ex_7_11.sce new file mode 100755 index 000000000..2e785541b --- /dev/null +++ b/50/CH7/EX7.11/ex_7_11.sce @@ -0,0 +1,32 @@ +// example 7.11 +// solve the boundary value problem u''=u'+1; +// u(0)=1; u(x=1)=2(%e-1); h=1/3; + + +// we know; u''=(u(j-1)-2*u(j)+u(j+1))/h^2; +// we know; u'=(u(j+1)-u(j-1))/2h; + +// 1) second order method; + x=0:1/3:1; + + u=[u0 u1 u2 u3 ]; +// hence; +disp('(u(j-1)-2*u(j)+u(j+1))/h^2=((u(j+1)-u(j-1))/2h)+1') // for j=1,2; + + +disp('for j=1 (7/6)*u0-2*u1+(5/6)*u2=(1/9)') + +disp('for j=2 (7/6)*u1-2*u2+(5/6)*u3=(1/9)') + + +// hence eliminating u1! +// solving for u1,u2, +u0=1; +u3=2*(%e-1); +u1=1.454869; +u2=2.225019; + +disp(x); +disp(u); + + diff --git a/50/CH7/EX7.3/ex_7_3.sce b/50/CH7/EX7.3/ex_7_3.sce new file mode 100755 index 000000000..6c8f5c1c5 --- /dev/null +++ b/50/CH7/EX7.3/ex_7_3.sce @@ -0,0 +1,29 @@ +// example 7.3 +// solve by shooting method; + +// u''=2*u*u'; +// u(0)=0.5; u(1)=1; + +// let -> U1(x)=du/dx; +// U2(x)=d2u/dx2; + +// U(x)=[U1(x);U2(x)] + +// hence ; +// dU/dx=f(x,U); + +h=.25; + +ub=[.5,1]; + +up=[0:.1:1]; + +x=0:h:1; + +deff('[w]=f(x,U)','w=[U(2); 2*U(1)*U(2)]') + + + +[U] = shooting(ub,up,x,f); + +// the solution obtained would show the values of u in the first collumn and their corresponding derivatives in the second collumn ; \ No newline at end of file diff --git a/50/CH7/EX7.4/ex_7_4.sce b/50/CH7/EX7.4/ex_7_4.sce new file mode 100755 index 000000000..2536351fd --- /dev/null +++ b/50/CH7/EX7.4/ex_7_4.sce @@ -0,0 +1,28 @@ +// example 7.4 +// solve by shooting method; + +// u''=2*u*u'; +// u(0)=0.5; u(1)=1; + +// let -> U1(x)=du/dx; +// U2(x)=d2u/dx2; + +// U(x)=[U1(x);U2(x)] + +// hence ; +// dU/dx=f(x,U); + +h=.25; + +ub=[.5,1]; + +up=[0:.1:1]; + +x=0:h:1; + +deff('[w]=f(x,U)','w=[U(2); 2*U(1)*U(2)]') + + +[U] = shooting(ub,up,x,f); + +// the solution obtained would show the values of u in the first collumn and their corresponding derivatives in the second collumn ; \ No newline at end of file diff --git a/50/CH7/EX7.5/ex_7_5.sce b/50/CH7/EX7.5/ex_7_5.sce new file mode 100755 index 000000000..8d6eed3f7 --- /dev/null +++ b/50/CH7/EX7.5/ex_7_5.sce @@ -0,0 +1,52 @@ +// example 7.5 +// solve the boundary value problem u''=u+x; +// u(x=0)=u(0)=0; u(x=1)=u(4)=0; h=1/4; + + +// we know; u''=(u(j-1)-2*u(j)+u(j+1))/h^2; + +// 1) second order method; + x=0:1/4:1; + u0=0; + u4=0; + u=[u0 u1 u2 u3 u4]; +// hence; +disp('(u(j-1)-2*u(j)+u(j+1))/h^2=u(j)+x(j)') // for j=1,2,3; + +disp('for j=1 -16*u0+33*u1-16*u2=-.25') + +disp('for j=2 -16*u1+33*u2-16*u3=-.50') + +disp('for j=3 -16*u2+33*u3-16*u4=-.75') + +// hence solving for u1,u2,u3) , +u1=-.034885; +u2=-.056326; +u3=-.050037; + +disp(x); +disp(u); + +// 2) numerov method; + x=0:1/4:1; + u0=0; + u4=0; + u=[u0 u1 u2 u3 u4]; +// since according to numerov method we get the following system of equations; +disp('(191*u(j-1)-394*u(j)+191*u(j+1)=x(j-1)+10*x(j)+x(j+1)') // for j=1,2,3; + +disp('for j=1 191*u0-394*u1+191*u2=3') + +disp('for j=2 191*u1-394*u2+191*u3=6') + +disp('for j=3 191*u2-394*u3+191*u4=9') + +// hence solving for u1,u2,u3 , +u1=-.034885 +u2=-.056326 +u3=-.050037 + + +disp(x); +disp(u); + diff --git a/50/CH7/EX7.6/ex_7_6.sce b/50/CH7/EX7.6/ex_7_6.sce new file mode 100755 index 000000000..04a28e398 --- /dev/null +++ b/50/CH7/EX7.6/ex_7_6.sce @@ -0,0 +1,32 @@ +// example 7.6 +// solve the boundary value problem u''=u*x; +// u(0)+u'(0)=1; u(x=1)=0; h=1/3; + + +// we know; u''=(u(j-1)-2*u(j)+u(j+1))/h^2; + +// 1) second order method; + x=0:1/3:1; + + u3=1; + u=[u0 u1 u2 u3 ]; +// hence; +disp('(u(j-1)-2*u(j)+u(j+1))/h^2=u(j)*x(j)') // for j=0,1,2,3; + +disp('for j=0 u1!-2*u0+u1=0') // u1!=u(-1) + +disp('for j=1 u0-2*u1+u2=(1/27)u1') + +disp('for j=2 u1-2*u2+u3=(2/27)u2') + +// we know; u'=(u(j+1)-u(j-1))/2h +// hence eliminating u1! +// solving for u0,u1,u2,u3 , +u0=-.9879518; +u1=-.3253012; +u2=-.3253012; + +disp(x); +disp(u); + + diff --git a/50/DEPENDENCIES/Euler1.sce b/50/DEPENDENCIES/Euler1.sce new file mode 100755 index 000000000..67fa05b6d --- /dev/null +++ b/50/DEPENDENCIES/Euler1.sce @@ -0,0 +1,26 @@ +function [u,t] = Euler1(u0,t0,tn,h,f) + +//Euler 1st order method solving ODE +// du/dt = f(u,t), with initial +//conditions u=u0 at t=t0. The +//solution is obtained for t = [t0:h:tn] +//and returned in u + +umaxAllowed = 1e+100; + +t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0; + +for j = 1:n-1 + u(j+1) = u(j) + h*f(t(j),u(j)); + if u(j+1) > umaxAllowed then + disp('Euler 1 - WARNING: underflow or overflow'); + disp('Solution sought in the following range:'); + disp([t0 h tn]); + disp('Solution evaluated in the following range:'); + disp([t0 h t(j)]); + n = j; t = t(1,1:n); u = u(1,1:n); + break; + end; +end; + +endfunction diff --git a/50/DEPENDENCIES/F28.sce b/50/DEPENDENCIES/F28.sce new file mode 100755 index 000000000..e79b88348 --- /dev/null +++ b/50/DEPENDENCIES/F28.sce @@ -0,0 +1,8 @@ +function [F] = f28(x) + +//f2(x) = 0, with x = [x(1);x(2)] +//represents a system of 2 non-linear equations +F(1) = x(1)^2+x(1)*x(2) + x(2)^2 - 7; +F(2) = x(1)^3+x(2)^3 -9; + +endfunction diff --git a/50/DEPENDENCIES/LandU.sce b/50/DEPENDENCIES/LandU.sce new file mode 100755 index 000000000..fee5d09ea --- /dev/null +++ b/50/DEPENDENCIES/LandU.sce @@ -0,0 +1,12 @@ +function [U,L]=LandU(A,n) + U=A + L=eye(n,n) + for p=1:1:n-1 + for i=p+1:1:n + m=A(i,p)/A(p,p); + L(i,p)=m; + A(i,:)=A(i,:)-m*A(p,:); + U=A; + end + end +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/NBDP.sci b/50/DEPENDENCIES/NBDP.sci new file mode 100755 index 000000000..535583f62 --- /dev/null +++ b/50/DEPENDENCIES/NBDP.sci @@ -0,0 +1,47 @@ +function [P]=NBDP(x,n,xL,f) +//This function calculates a Newton Forward-Difference Polynomial of +//order n, evaluated at x, using column vectors xL, f as the reference +//table. The first value of xL and of f, represent, respectively, +//xo and fo in the equation for the polynomial. +[m,nc]=size(f) +//check that it is indeed a column vector +if (nc<>1) then + error('f is not a column vector.'); + abort +end; +//check the difference order +if (n >= m) then + disp(n,"n="); + disp(m,"m="); + error('n must be less than or equal to m-1'); + abort +end; +// +xo = xL(m,1); +delx = mtlb_diff(xL); +h = delx(1,1); +s = (x-xo)/h; +P = f(m,1); +delf = f; +disp(delf); +for i = 1:n + delf = mtlb_diff(delf); + [m,nc] = size(delf); + disp(delf); + P = P + Binomial(s+i-1,i)*delf(m,1) +end; +endfunction + +function[C]=Binomial(s,i) + C = 1.0; + for k = 0:i-1 + C = C*(s-k); + end; + C = C/factorial(i) +endfunction +function[fact]=factorial(nn) + fact = 1.0 + for k = nn:-1:1 + fact=fact*k + end; +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/NDDinterpol.sci b/50/DEPENDENCIES/NDDinterpol.sci new file mode 100755 index 000000000..5b49ebc9c --- /dev/null +++ b/50/DEPENDENCIES/NDDinterpol.sci @@ -0,0 +1,17 @@ +function P1=NDDinterpol (x0,x1,f0,f1) + x=poly(0,"x"); + f01=(f1-f0)/(x1-x0); + P1=f0+(x-x0)*f01; +endfunction + + + + + + + + + + + + \ No newline at end of file diff --git a/50/DEPENDENCIES/NDDinterpol2.sci b/50/DEPENDENCIES/NDDinterpol2.sci new file mode 100755 index 000000000..9618c8e9e --- /dev/null +++ b/50/DEPENDENCIES/NDDinterpol2.sci @@ -0,0 +1,19 @@ +function P2=NDDinterpol2 (x,f) + X=poly(0,"X"); + f01=(f(2)-f(1))/(x(2)-x(1)); + f13=(f(3)-f(2))/(x(3)-x(2)); + f013=(f13-f01)/(x(3)-x(1)); + P2=f(1)+(X-x(1))*f01+(X-x(1))*(X-x(2))*f013; +endfunction + + + + + + + + + + + + \ No newline at end of file diff --git a/50/DEPENDENCIES/RK4.sci b/50/DEPENDENCIES/RK4.sci new file mode 100755 index 000000000..7539a4b15 --- /dev/null +++ b/50/DEPENDENCIES/RK4.sci @@ -0,0 +1,30 @@ +function [u,t] = RK4(u0,t0,tn,h,f) + +// RK4 method solving ODE +// du/dt = f(u,t), with initial +//conditions u=u0 at t=t0. The +//solution is obtained for t = [t0:h:tn] +//and returned in u + +umaxAllowed = 1e+100; + +t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0; + +for j = 1:n-1 + k1=h*f(t(j),u(j)); + k2=h*f(t(j)+h/2,u(j)+k1/2); + k3=h*f(t(j)+h/2,u(j)+k2/2); + k4=h*f(t(j)+h,u(j)+k3); + u(j+1) = u(j) + (1/6)*(k1+2*k2+2*k3+k4); + if u(j+1) > umaxAllowed then + disp('Euler 1 - WARNING: underflow or overflow'); + disp('Solution sought in the following range:'); + disp([t0 h tn]); + disp('Solution evaluated in the following range:'); + disp([t0 h t(j)]); + n = j; t = t(1,1:n); u = u(1,1:n); + break; + end; +end; + +endfunction diff --git a/50/DEPENDENCIES/Vaitken.sce b/50/DEPENDENCIES/Vaitken.sce new file mode 100755 index 000000000..cf93e1b3a --- /dev/null +++ b/50/DEPENDENCIES/Vaitken.sce @@ -0,0 +1,9 @@ +// this program is exclusively coded to perform one iteration of aitken method, + +function x0aa=aitken(x0,x1,x2,g) +x0a=x0-(x1-x0)^2/(x2-2*x1+x0); +x1a=g(x0a); +x2a=g(x1a); +x0aa=x0a-(x1a-x0a)^2/(x2a-2*x1a+x0a); + +endfunction diff --git a/50/DEPENDENCIES/Vbisection.sce b/50/DEPENDENCIES/Vbisection.sce new file mode 100755 index 000000000..cb23696bf --- /dev/null +++ b/50/DEPENDENCIES/Vbisection.sce @@ -0,0 +1,28 @@ +function x=bisection(a,b,f) + N=100; //define max. number of iterations + PE=10^-4 //define tolerance + if (f(a)*f(b) > 0) then + error('no root possible f(a)*f(b) > 0') // checking if the decided range is containing a root + abort; + end; + if(abs(f(a)) maxval) then error('Solution diverges'); + abort; + break + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/Vgausselim.sce b/50/DEPENDENCIES/Vgausselim.sce new file mode 100755 index 000000000..9488608a1 --- /dev/null +++ b/50/DEPENDENCIES/Vgausselim.sce @@ -0,0 +1,44 @@ +function [x] = gausselim(A,b) + +//This function obtains the solution to the system of +//linear equations A*x = b, given the matrix of coefficients A +//and the right-hand side vector, b + +[nA,mA] = size(A) +[nb,mb] = size(b) + +if nA<>mA then + error('gausselim - Matrix A must be square'); + abort; +elseif mA<>nb then + error('gausselim - incompatible dimensions between A and b'); + abort; +end; + +a = [A b]; + +//Forward elimination + +n = nA; +for k=1:n-1 + for i=k+1:n + for j=k+1:n+1 + a(i,j)=a(i,j)-a(k,j)*a(i,k)/a(k,k); + end; + end; +end; + +//Backward substitution + +x(n) = a(n,n+1)/a(n,n); + +for i = n-1:-1:1 + sumk=0 + for k=i+1:n + sumk=sumk+a(i,k)*x(k); + end; + x(i)=(a(i,n+1)-sumk)/a(i,i); +end; + +endfunction + diff --git a/50/DEPENDENCIES/Vgaussseidel.sce b/50/DEPENDENCIES/Vgaussseidel.sce new file mode 100755 index 000000000..ba227af87 --- /dev/null +++ b/50/DEPENDENCIES/Vgaussseidel.sce @@ -0,0 +1,24 @@ +function [X]=gaussseidel(A,n,N,X,b) + L=A; + U=A; + D=A; + for i=1:1:n + for j=1:1:n + if j>i then L(i,j)=0; + D(i,j)=0; + end + if i>j then U(i,j)=0; + D(i,j)=0; + end + if i==j then L(i,j)=0; + U(i,j)=0; + end + end + + end + for k=1:1:N + X=(D+L)^-1*(-U*X+b); + disp(X) + end + +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/Vnewton.sce b/50/DEPENDENCIES/Vnewton.sce new file mode 100755 index 000000000..f0fb337be --- /dev/null +++ b/50/DEPENDENCIES/Vnewton.sce @@ -0,0 +1,17 @@ + +function x=newton(x,f,fp) + R=100; + PE=10^-8; + maxval=10^4; + + for n=1:1:R + x=x-f(x)/fp(x); + if abs(f(x))<=PE then break + end + if (abs(f(x))>maxval) then error('Solution diverges'); + abort + break + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/Vnewton4.sce b/50/DEPENDENCIES/Vnewton4.sce new file mode 100755 index 000000000..d35f210b4 --- /dev/null +++ b/50/DEPENDENCIES/Vnewton4.sce @@ -0,0 +1,17 @@ +function x=newton4(x,f,fp) + R=4; + PE=10^-15; + maxval=10^4; + for n=1:1:R + if fp(x)==0 then disp("select another initial root x0") + end + x=x-f(x)/fp(x); + if abs(f(x))<=PE then break + end + if (abs(f(x))>maxval) then error('Solution diverges'); + abort + break + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/Vregulafalsi.sce b/50/DEPENDENCIES/Vregulafalsi.sce new file mode 100755 index 000000000..7b6b0259a --- /dev/null +++ b/50/DEPENDENCIES/Vregulafalsi.sce @@ -0,0 +1,12 @@ +function [x]=regularfalsi(a,b,f) + N=100; + PE=10^-5; + for n=2:1:N + x=a-(a-b)*f(a)/(f(a)-f(b)); + if abs(f(x))<=PE then break; + elseif (f(a)*f(x)<0) then b=x; + else a=x; + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/Vsecant.sce b/50/DEPENDENCIES/Vsecant.sce new file mode 100755 index 000000000..a44c89157 --- /dev/null +++ b/50/DEPENDENCIES/Vsecant.sce @@ -0,0 +1,12 @@ +function [x]=secant(a,b,f) + N=100; // define max. no. iterations to be performed + PE=10^-4 // define tolerance for convergence + for n=1:1:N // initiating for loop + x=a-(a-b)*f(a)/(f(a)-f(b)); + if abs(f(x))<=PE then break; //checking for the required condition + else a=b; + b=x; + end + end + disp(n," no. of iterations =") // +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/Vsim_eulercauchy.sce b/50/DEPENDENCIES/Vsim_eulercauchy.sce new file mode 100755 index 000000000..129472469 --- /dev/null +++ b/50/DEPENDENCIES/Vsim_eulercauchy.sce @@ -0,0 +1,24 @@ +function [u,v,t] = simeulercauchy(u0,v0,t0,tn,h,f1,f2) + + +// du/dt = f1(t,u,v), dv/dt = f2(t,u,v) with initial +//conditions u=u0,v=v0 at t=t0. The +//solution is obtained for t = [t0:h:tn] +//and returned in u,v + + +umaxAllowed = 1e+100; + +t = [t0:h:tn]; u = zeros(t);v= zeros(t); n = length(u); u(1) = u0;v(1)=v0; + +for j = 1:n-1 + k11=h*f1(t(j),u(j),v(j)); + k21=h*f2(t(j),u(j),v(j)); + k12=h*f1(t(j)+h,u(j)+k11,v(j)+k21); + k22=h*f2(t(j)+h,u(j)+k11,v(j)+k21); + u(j+1) = u(j) + (k11+k12)/2; + v(j+1) = v(j) + (k21+k22)/2; + +end; + +endfunction diff --git a/50/DEPENDENCIES/adamsbashforth3.sci b/50/DEPENDENCIES/adamsbashforth3.sci new file mode 100755 index 000000000..dfc65961f --- /dev/null +++ b/50/DEPENDENCIES/adamsbashforth3.sci @@ -0,0 +1,23 @@ +function [u,t] = adamsbashforth3(u0,t0,tn,h,f) + +//adamsbashforth3 3rd order method solving ODE +// du/dt = f(u,t), with initial +//conditions u=u0 at t=t0. The +//solution is obtained for t = [t0:h:tn] +//and returned in u + +umaxAllowed = 1e+100; + +t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0; +for j = 1:n-2 +if j<3 then + k1=h*f(t(j),u(j)); + k2=h*f(t(j)+h,u(j)+k1); + u(j+1) = u(j) + (k2+k1)/2; +end; + +if j>=2 then + u(j+2) = u(j+1) + (h/12)*(23*f(t(j+1),u(j+1))-16*f(t(j),u(j))+5*f(t(j-1),u(j-1))); +end; +end; +endfunction diff --git a/50/DEPENDENCIES/aitken.sce b/50/DEPENDENCIES/aitken.sce new file mode 100755 index 000000000..cf93e1b3a --- /dev/null +++ b/50/DEPENDENCIES/aitken.sce @@ -0,0 +1,9 @@ +// this program is exclusively coded to perform one iteration of aitken method, + +function x0aa=aitken(x0,x1,x2,g) +x0a=x0-(x1-x0)^2/(x2-2*x1+x0); +x1a=g(x0a); +x2a=g(x1a); +x0aa=x0a-(x1a-x0a)^2/(x2a-2*x1a+x0a); + +endfunction diff --git a/50/DEPENDENCIES/aitkeninterpol.sci b/50/DEPENDENCIES/aitkeninterpol.sci new file mode 100755 index 000000000..2b6052d81 --- /dev/null +++ b/50/DEPENDENCIES/aitkeninterpol.sci @@ -0,0 +1,16 @@ +function P1=aitkeninterpol (x0,x1,f0,f1) + x=poly(0,"x"); + P1=(1/(x1-x0))*det([f0 x0-x;f1 x1-x]); +endfunction + + + + + + + + + + + + \ No newline at end of file diff --git a/50/DEPENDENCIES/back.sce b/50/DEPENDENCIES/back.sce new file mode 100755 index 000000000..d92199af3 --- /dev/null +++ b/50/DEPENDENCIES/back.sce @@ -0,0 +1,16 @@ +function [x] = back(U,Z) + +x=zeros(1,n); +for i = n:-1:1 + sumk=0 + for j=i+1:n + sumk=sumk+U(i,j)*x(j); + end; + x(i)=(Z(i)-sumk)/U(i,i); +end; + + + + +endfunction + diff --git a/50/DEPENDENCIES/backeuler.sci b/50/DEPENDENCIES/backeuler.sci new file mode 100755 index 000000000..714f130ce --- /dev/null +++ b/50/DEPENDENCIES/backeuler.sci @@ -0,0 +1,22 @@ +function [u] = backeuler(u0,t0,tn,h,f) + +//backeuler 1st order method solving ODE +// du/dt = f(u,t), with initial +//conditions u=u0 at t=t0. The +//solution is obtained for t = [t0:h:tn] +//and returned in u + +umaxAllowed = 1e+100; + + t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0; + +for j=1:n-1 + u(j+1)=u(j); +for i = 0:5 + u(j+1) = u(j) + h*f(t(j+1),u(j+1)); + i=i+1; +end; +end; + + +endfunction diff --git a/50/DEPENDENCIES/chebyshev.sce b/50/DEPENDENCIES/chebyshev.sce new file mode 100755 index 000000000..b883af5fa --- /dev/null +++ b/50/DEPENDENCIES/chebyshev.sce @@ -0,0 +1,18 @@ +function x=chebyshev(x,f,fp,fpp) + R=100; + PE=10^-5; + maxval=10^4; + if fp(x)==0 then disp("select another initial root x0"); + break; + end + for n=1:1:R + x=x-f(x)/fp(x)-(1/2)*(f(x)/fp(x))^2 *(fpp(x)/fp(x)); + if abs(f(x))<=PE then break; + end + if (abs(f(x))>maxval) then error('Solution diverges'); + abort; + break + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/cholesky.sce b/50/DEPENDENCIES/cholesky.sce new file mode 100755 index 000000000..b3911d495 --- /dev/null +++ b/50/DEPENDENCIES/cholesky.sce @@ -0,0 +1,16 @@ +function L=cholesky (A,n) + L=zeros(n,n); + for k=1:1:n + S=0; + P=0; + for j=1:1:k-1 + S=S+(L(k,j)^2); + P=P+L(i,j)*L(k,j) + end + L(k,k)=sqrt(A(k,k)-S); + for i=k+1:1:n + L(i,k)=(A(i,k)-P)/L(k,k); + end + end + +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/comp_trapezoidal.sci b/50/DEPENDENCIES/comp_trapezoidal.sci new file mode 100755 index 000000000..f76a78571 --- /dev/null +++ b/50/DEPENDENCIES/comp_trapezoidal.sci @@ -0,0 +1,34 @@ +function I=comptrapezoidal(x,h,f) + //This function calculates the numerical integration of f(x)dx +//between limits x(1) and x(n) using composite trapezoidal rule +//Check that x and y have the same size (which must be an odd number) +//Also, the values of x must be equally spaced with spacing h +y=feval(x,f); +[nrx,ncx]=size(x) +[nrf,ncf]=size(y) +if ((nrx<>1)|(nrf<>1)) then + error('x or f, or both, not column vector(s)'); + abort; +end; +if ((ncx<>ncf)) then + error('x and f are not of the same length'); + abort; +end; +//check that the size of the lists xL and f is odd +if (modulo(ncx,2)==0) then + disp(ncx,"list size =") + error('list size must be an odd number'); + abort +end; +n = ncx; + +I = f(x(1)) + f(x(n)); +for j = 2:n-1 + if(modulo(j,2)==0) then + I = I + 2*f(x(j)); + else + I = I + 2*f(x(j)); + end; +end; +I = (h/2.0)*I +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/eigenvectors.sci b/50/DEPENDENCIES/eigenvectors.sci new file mode 100755 index 000000000..e36ffaaeb --- /dev/null +++ b/50/DEPENDENCIES/eigenvectors.sci @@ -0,0 +1,37 @@ +function [x,lam] = geigenvectors(A,B) + +//Calculates unit eigenvectors of matrix A +//returning a matrix x whose columns are +//the eigenvectors. The function also +//returns the eigenvalues of the matrix. + +[nA,mA] = size(A); +[nB,mB] = size(B); + +if (mA<>nA | mB<>nB) then + error('geigenvectors - matrix A or B not square'); + abort; +end; + +if nA<>nB then + error('geigenvectors - matrix A and B have different dimensions'); + abort; +end; + +lam = poly(0,'lam'); //Define variable "lam" +chPoly = det(A-B*lam); //Characteristic polynomial +lam = roots(chPoly)'; //Eigenvalues of matrix A + +x = []; n = nA; + +for k = 1:n + BB = A - lam(k)*B; //Characteristic matrix + CC = BB(1:n-1,1:n-1); //Coeff. matrix for reduced system + bb = -BB(1:n-1,n); //RHS vector for reduced system + y = CC\bb; //Solution for reduced system + y = [y;1]; //Complete eigenvector + + x = [x y]; //Add eigenvector to matrix +end; + +endfunction diff --git a/50/DEPENDENCIES/eulermidpoint.sci b/50/DEPENDENCIES/eulermidpoint.sci new file mode 100755 index 000000000..f05453299 --- /dev/null +++ b/50/DEPENDENCIES/eulermidpoint.sci @@ -0,0 +1,26 @@ +function [u] = eulermidpoint(u0,t0,tn,h,f,fp) + +//midpoint 1st order method solving ODE +// du/dt = f(u,t), with initial +//conditions u=u0 at t=t0. The +//solution is obtained for t = [t0:h:tn] +//and returned in u + +umaxAllowed = 1e+100; + +t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0; +u(2)=u(1)+h*f(t(1),u(1))+(h^2/2)*fp(t(1),u(1)); +for j = 2:n-1 + u(j+1) = u(j-1) + 2*h*f(t(j),u(j)); + if u(j+1) > umaxAllowed then + disp('Euler 1 - WARNING: underflow or overflow'); + disp('Solution sought in the following range:'); + disp([t0 h tn]); + disp('Solution evaluated in the following range:'); + disp([t0 h t(j)]); + n = j; t = t(1,1:n); u = u(1,1:n); + break; + end; +end; + +endfunction diff --git a/50/DEPENDENCIES/fact.sci b/50/DEPENDENCIES/fact.sci new file mode 100755 index 000000000..c12441e95 --- /dev/null +++ b/50/DEPENDENCIES/fact.sci @@ -0,0 +1,6 @@ +function x=fact(n) + x=1; + for i=2:1:n + x=x*i; + end; +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/fore.sce b/50/DEPENDENCIES/fore.sce new file mode 100755 index 000000000..d8338b7cb --- /dev/null +++ b/50/DEPENDENCIES/fore.sce @@ -0,0 +1,12 @@ +function x=fore(L,b) + +for i = 1:1:n + sumk=0 + for j=1:i-1 + sumk=sumk+L(i,j)*x(j); + end; + x(i)=(b(i)-sumk)/L(i,i); +end; + +endfunction + diff --git a/50/DEPENDENCIES/geigenvectors.sci b/50/DEPENDENCIES/geigenvectors.sci new file mode 100755 index 000000000..e36ffaaeb --- /dev/null +++ b/50/DEPENDENCIES/geigenvectors.sci @@ -0,0 +1,37 @@ +function [x,lam] = geigenvectors(A,B) + +//Calculates unit eigenvectors of matrix A +//returning a matrix x whose columns are +//the eigenvectors. The function also +//returns the eigenvalues of the matrix. + +[nA,mA] = size(A); +[nB,mB] = size(B); + +if (mA<>nA | mB<>nB) then + error('geigenvectors - matrix A or B not square'); + abort; +end; + +if nA<>nB then + error('geigenvectors - matrix A and B have different dimensions'); + abort; +end; + +lam = poly(0,'lam'); //Define variable "lam" +chPoly = det(A-B*lam); //Characteristic polynomial +lam = roots(chPoly)'; //Eigenvalues of matrix A + +x = []; n = nA; + +for k = 1:n + BB = A - lam(k)*B; //Characteristic matrix + CC = BB(1:n-1,1:n-1); //Coeff. matrix for reduced system + bb = -BB(1:n-1,n); //RHS vector for reduced system + y = CC\bb; //Solution for reduced system + y = [y;1]; //Complete eigenvector + + x = [x y]; //Add eigenvector to matrix +end; + +endfunction diff --git a/50/DEPENDENCIES/generaliteration.sce b/50/DEPENDENCIES/generaliteration.sce new file mode 100755 index 000000000..3e89287e1 --- /dev/null +++ b/50/DEPENDENCIES/generaliteration.sce @@ -0,0 +1,21 @@ + +function x=generaliteration(x,g,gp) + R=5; + PE=10^-8; + maxval=10^4; + k=gp(x); + if abs(k)>1 then error('function chosen does not converge') + abort; + end + for n=1:1:R + x=g(x); + disp(x); + if abs(g(x))<=PE then break + end + if (abs(g(x))>maxval) then error('Solution diverges'); + abort + break + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/generaliteration2.sce b/50/DEPENDENCIES/generaliteration2.sce new file mode 100755 index 000000000..29545493a --- /dev/null +++ b/50/DEPENDENCIES/generaliteration2.sce @@ -0,0 +1,22 @@ + +function x=generaliteration2(x,g,gp) + R=2; + PE=10^-8; + maxval=10^4; + A=[0 0]; + k=gp(x); + if abs(k)>1 then error('function chosen does not converge') + abort; + end + for n=1:1:R + x=g(x); + disp(x); + if abs(g(x))<=PE then break + end + if (abs(g(x))>maxval) then error('Solution diverges'); + abort + break + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/hermiteinterpol.sci b/50/DEPENDENCIES/hermiteinterpol.sci new file mode 100755 index 000000000..d9bc22148 --- /dev/null +++ b/50/DEPENDENCIES/hermiteinterpol.sci @@ -0,0 +1,37 @@ +function P= hermiteinterpol(x,f,fp) + X=poly(0,"X"); + function f0=L0(X) + + f0=(X-x(2))*(X-x(3))/((x(1)-x(2))*(x(1)-x(3))) + endfunction; + a0=[1-2*(X-x(1))*numdiff(L0,x(1))]; +L0=(X-x(2))*(X-x(3))/((x(1)-x(2))*(x(1)-x(3))); + A0=a0*L0*L0; + disp(A0) + B0=(X-x(1))*L0^2; + + X=poly(0,"X"); + + function f1=L1(X) + + f1=(X-x(1))*(X-x(3))/((x(2)-x(1))*(x(2)-x(3))) + endfunction; +a1=[1-2*(X-x(2))*0]; +L1=(X-x(1))*(X-x(3))/((x(2)-x(1))*(x(2)-x(3))); +A1=a1 *L1*L1; +disp(A1) + B1=(X-x(2))*L1^2; + function f2=L2(X) + + f2=(X-x(1))*(X-x(2))/((x(3)-x(1))*(x(3)-x(2))) + endfunction; +a2=[1-2*(X-x(3))*numdiff(L2,x(3))]; +L2=(X-x(1))*(X-x(2))/((x(3)-x(1))*(x(3)-x(2))); +A2=a2 *L2*L2; +disp(A2) + B2=(X-x(3))*L2^2; + + + + P=A0*f(1)+A1*f(2)+A2*f(3)+B0*fp(1)+B1*fp(2)+B2*fp(3); +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/heun.sci b/50/DEPENDENCIES/heun.sci new file mode 100755 index 000000000..f308280b0 --- /dev/null +++ b/50/DEPENDENCIES/heun.sci @@ -0,0 +1,28 @@ +function [u,t] = heun(u0,t0,tn,h,f) + +//heun method solving ODE +// du/dt = f(u,t), with initial +//conditions u=u0 at t=t0. The +//solution is obtained for t = [t0:h:tn] +//and returned in u + +umaxAllowed = 1e+100; + +t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0; + +for j = 1:n-1 + k1=h*f(t(j),u(j)); + k2=h*f(t(j)+h,u(j)+k1); + u(j+1) = u(j) + (k2+k1)/2; + if u(j+1) > umaxAllowed then + disp('Euler 1 - WARNING: underflow or overflow'); + disp('Solution sought in the following range:'); + disp([t0 h tn]); + disp('Solution evaluated in the following range:'); + disp([t0 h t(j)]); + n = j; t = t(1,1:n); u = u(1,1:n); + break; + end; +end; + +endfunction diff --git a/50/DEPENDENCIES/iteratedinterpol.sci b/50/DEPENDENCIES/iteratedinterpol.sci new file mode 100755 index 000000000..2997bd6c5 --- /dev/null +++ b/50/DEPENDENCIES/iteratedinterpol.sci @@ -0,0 +1,18 @@ +function [L012,L02,L01]=iteratedinterpol (x,f,n) + X=poly(0,"X"); + L01=(1/(x(2)-x(1)))*det([f(1) x(1)-X;f(2) x(2)-X]); + L02=(1/(x(3)-x(1)))*det([f(1) x(1)-X;f(3) x(3)-X]); + L012=(1/(x(3)-x(2)))*det([L01 x(2)-X;L02 x(3)-X]); + +endfunction + + + + + + + + + + + \ No newline at end of file diff --git a/50/DEPENDENCIES/jacob28.sce b/50/DEPENDENCIES/jacob28.sce new file mode 100755 index 000000000..0553cd801 --- /dev/null +++ b/50/DEPENDENCIES/jacob28.sce @@ -0,0 +1,9 @@ +function [J] = jacob28(x) + +//Evaluates the Jacobian of a 2x2 +//system of non-linear equations + +J(1,1) = 2*x(1)+x(2); J(1,2) = x(1)+2*x(2); +J(2,1) = 3*x(1)^2; J(2,2) =3* x(2)^2; +endfunction + diff --git a/50/DEPENDENCIES/jacobianmat.sci b/50/DEPENDENCIES/jacobianmat.sci new file mode 100755 index 000000000..1a6d8ba75 --- /dev/null +++ b/50/DEPENDENCIES/jacobianmat.sci @@ -0,0 +1,8 @@ +function J= jacobianmat (f1,f2,h,k) + J=zeros(2,2); +J(1,1)=(f1(1+h,1)-f1(1,1))/2*h; + +J(1,2)=(f1(1,1+k)-f1(1,1))/2*k; +J(2,1)=(f2(1+h,1)-f2(1,1))/2*h; +J(2,2)=(f2(1,1+k)-f2(1,1))/2*k; +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/jacobiiteration.sce b/50/DEPENDENCIES/jacobiiteration.sce new file mode 100755 index 000000000..2cdc4e224 --- /dev/null +++ b/50/DEPENDENCIES/jacobiiteration.sce @@ -0,0 +1,23 @@ +function [X]=jacobiiteration(A,n,N,X,b) + L=A; + U=A; + D=A; + for i=1:1:n + for j=1:1:n + if j>i then L(i,j)=0; + D(i,j)=0; + end + if i>j then U(i,j)=0; + D(i,j)=0; + end + if i==j then L(i,j)=0; + U(i,j)=0; + end + end + + end + for k=1:1:N + X=-D^-1*(L+U)*X+D^-1*(b); + end + +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/jordan.sce b/50/DEPENDENCIES/jordan.sce new file mode 100755 index 000000000..22ba3c1fa --- /dev/null +++ b/50/DEPENDENCIES/jordan.sce @@ -0,0 +1,20 @@ +function [M]=jorden(A,b) + M=[A b]; + [ra,ca]=size(A); + [rb,cb]=size(b); + n=ra; + for p=1:1:n + for k=(p+1):1:n + if abs(M(k,p))>abs(M(p,p)) then + M({p,k},:)=M({k,p},:); + end + end + M(p,:)=M(p,:)/M(p,p); + for i=1:1:p-1 + M(i,:)=M(i,:)-M(p,:)*(M(i,p)/M(p,p)); + end + for i=p+1:1:n + M(i,:)=M(i,:)-M(p,:)*(M(i,p)/M(p,p)); + end + end +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/lagrangefundamentalpoly.sci b/50/DEPENDENCIES/lagrangefundamentalpoly.sci new file mode 100755 index 000000000..b55a42ff6 --- /dev/null +++ b/50/DEPENDENCIES/lagrangefundamentalpoly.sci @@ -0,0 +1,30 @@ +function P2=lagrangefundamentalpoly(x,f,n) + [nrx,ncx]=size(x) + [nrf,ncf]=size(f) +if ((nrx<>1)|(nrf<>1)) then + error('x or f, or both, not column vector(s)'); + abort; +end; +if ((ncx<>ncf)) then + error('x and f are not of the same length'); + abort; +end; + + X=poly(0,"X"); +L=zeros(n); + +P2=0; +for i=1:n+1 + L(i)=1; + for j=1:n+1 + if i~=j then + L(i)=L(i)*(X-x(j))/(x(i)-x(j)) + end; +end; +P2=P2+L(i)*f(i); +end; + + +endfunction + + diff --git a/50/DEPENDENCIES/legrangeinterpol.sci b/50/DEPENDENCIES/legrangeinterpol.sci new file mode 100755 index 000000000..cbaff9113 --- /dev/null +++ b/50/DEPENDENCIES/legrangeinterpol.sci @@ -0,0 +1,18 @@ +function P1=legrangeinterpol (x0,x1,f0,f1) + x=poly(0,"x"); + L0=(x-x1)/(x0-x1); + L1=(x-x0)/(x1-x0); + P1=L0*f0+L1*f1; +endfunction + + + + + + + + + + + + \ No newline at end of file diff --git a/50/DEPENDENCIES/linearinterpol.sci b/50/DEPENDENCIES/linearinterpol.sci new file mode 100755 index 000000000..a6c10cbc5 --- /dev/null +++ b/50/DEPENDENCIES/linearinterpol.sci @@ -0,0 +1,3 @@ +function fp=linearinterpol(xL,f) + fp=(f(2)-f(1))/(xL(2)-xL(1)); +endfunction; \ No newline at end of file diff --git a/50/DEPENDENCIES/modified_newton.sce b/50/DEPENDENCIES/modified_newton.sce new file mode 100755 index 000000000..29d60f47b --- /dev/null +++ b/50/DEPENDENCIES/modified_newton.sce @@ -0,0 +1,17 @@ + +function x=modified_newton(x,f,fp) + R=100; + PE=10^-8; + maxval=10^4; + + for n=1:1:R + x=x-m*f(x)/fp(x); + if abs(f(x))<=PE then break + end + if (abs(f(x))>maxval) then error('Solution diverges'); + abort + break + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/modifiedeuler.sci b/50/DEPENDENCIES/modifiedeuler.sci new file mode 100755 index 000000000..a515eff0e --- /dev/null +++ b/50/DEPENDENCIES/modifiedeuler.sci @@ -0,0 +1,28 @@ +function [u,t] = modifiedeuler(u0,t0,tn,h,f) + +//modifiedeuler 1st order method solving ODE +// du/dt = f(u,t), with initial +//conditions u=u0 at t=t0. The +//solution is obtained for t = [t0:h:tn] +//and returned in u + +umaxAllowed = 1e+100; + +t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0; + +for j = 1:n-1 + k1=h*f(t(j),u(j)); + k2=h*f(t(j)+h/2,u(j)+k1/2); + u(j+1) = u(j) + k2; + if u(j+1) > umaxAllowed then + disp('Euler 1 - WARNING: underflow or overflow'); + disp('Solution sought in the following range:'); + disp([t0 h tn]); + disp('Solution evaluated in the following range:'); + disp([t0 h t(j)]); + n = j; t = t(1,1:n); u = u(1,1:n); + break; + end; +end; + +endfunction diff --git a/50/DEPENDENCIES/muller3.sce b/50/DEPENDENCIES/muller3.sce new file mode 100755 index 000000000..75007237c --- /dev/null +++ b/50/DEPENDENCIES/muller3.sce @@ -0,0 +1,32 @@ +function x=muller3(x0,x1,x2,f) + R=3; + PE=10^-8; + maxval=10^4; + for n=1:1:R + + La=(x2-x1)/(x1-x0); + Da=1+La; + ga=La^2*f(x0)-Da^2*f(x1)+(La+Da)*f(x2); + Ca=La*(La*f(x0)-Da*f(x1)+f(x2)); + + q=ga^2-4*Da*Ca*f(x2); + if q<0 then q=0; + end + p= sqrt(q); + if ga<0 then p=-p; + end + La=-2*Da*f(x2)/(ga+p); + x=x2+(x2-x1)*La; + if abs(f(x))<=PE then break + end + if (abs(f(x))>maxval) then error('Solution diverges'); + abort; + break + else + x0=x1; + x1=x2; + x2=x; + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/muller5.sce b/50/DEPENDENCIES/muller5.sce new file mode 100755 index 000000000..e2cf73b8a --- /dev/null +++ b/50/DEPENDENCIES/muller5.sce @@ -0,0 +1,32 @@ +function x=muller5(x0,x1,x2,f) + R=5; + PE=10^-8; + maxval=10^4; + for n=1:1:R + + La=(x2-x1)/(x1-x0); + Da=1+La; + ga=La^2*f(x0)-Da^2*f(x1)+(La+Da)*f(x2); + Ca=La*(La*f(x0)-Da*f(x1)+f(x2)); + + q=ga^2-4*Da*Ca*f(x2); + if q<0 then q=0; + end + p= sqrt(q); + if ga<0 then p=-p; + end + La=-2*Da*f(x2)/(ga+p); + x=x2+(x2-x1)*La; + if abs(f(x))<=PE then break + end + if (abs(f(x))>maxval) then error('Solution diverges'); + abort; + break + else + x0=x1; + x1=x2; + x2=x; + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/multipoint_iteration31.sce b/50/DEPENDENCIES/multipoint_iteration31.sce new file mode 100755 index 000000000..5c016939d --- /dev/null +++ b/50/DEPENDENCIES/multipoint_iteration31.sce @@ -0,0 +1,14 @@ +function x=multipoint_iteration31(x,f,fp,R) + R=3; + PE=10^-5; + maxval=10^4; + for n=1:1:R + x=x-f(x)/fp(x-(1/2)*(f(x)/fp(x))); + if abs(f(x))<=PE then break; + end + if (abs(f(x))>maxval) then error('Solution diverges'); + break + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/multipoint_iteration33.sce b/50/DEPENDENCIES/multipoint_iteration33.sce new file mode 100755 index 000000000..bd4bc3308 --- /dev/null +++ b/50/DEPENDENCIES/multipoint_iteration33.sce @@ -0,0 +1,14 @@ +function x=multipoint_iteration33(x,f,fp,R) + R=3; + PE=10^-5; + maxval=10^4; + for n=1:1:R + x=x-f(x)/fp(x)-f(x-(f(x)/fp(x)))/fp(x); + if abs(f(x))<=PE then break; + end + if (abs(f(x))>maxval) then error('Solution diverges'); + break + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/newton2.sce b/50/DEPENDENCIES/newton2.sce new file mode 100755 index 000000000..12714759c --- /dev/null +++ b/50/DEPENDENCIES/newton2.sce @@ -0,0 +1,21 @@ +function [x]=newton2(x0,f,fp) +//newton-raphson algorithm +N = 100; eps = 1.e-5; // define max. no. iterations and error +maxval = 10000.0; // define value for divergence +xx = x0; +while (N>0) + xn = xx-f(xx)/fp(xx); + if(abs(f(xn))<=eps) then x=xn; + disp(100-N); + return(x); + end; + if (abs(f(xx))>maxval) then disp(100-N); + error('Solution diverges'); + abort; + end; + N = N - 1; + xx = xn; +end; +error('No convergence'); +abort; +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/newton63.sce b/50/DEPENDENCIES/newton63.sce new file mode 100755 index 000000000..186634d2e --- /dev/null +++ b/50/DEPENDENCIES/newton63.sce @@ -0,0 +1,17 @@ + +function x=newton63(x,f,fp,fpp) + R=100; + PE=10^-15; + maxval=10^4; + + for n=1:1:R + x=x-(f(x)*fp(x))/(fp(x)^2-f(x)*fpp(x)); + if abs(f(x))<=PE then break + end + if (abs(f(x))>maxval) then error('Solution diverges'); + abort + break + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/newtonNLE.sce b/50/DEPENDENCIES/newtonNLE.sce new file mode 100755 index 000000000..cb0e9d5c9 --- /dev/null +++ b/50/DEPENDENCIES/newtonNLE.sce @@ -0,0 +1,31 @@ +function [xn] =newtonNLE(x0,f,J) + +//Newton-Raphson method applied to a +//system of linear equations f(x) = 0, +//given the jacobian function J, with +//J is the jacobian matrix of the functions. +//x = [x1;x2;...;xn], f = [f1;f2;...;fn] +//x0 is an initial guess of the solution + +N = 3; //define max. number of iterations +PE= 10^-15; //define tolerance +maxval = 10000.0; //define value for divergence +xx = x0; //load initial guess +for n=1:1:N + JJ = J(xx); + if(abs(det(JJ))=maxval) then + error('Solution diverges'); + abort; + end; + +end; +disp(n," no. of iterations ="); +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/newtonrap.sce b/50/DEPENDENCIES/newtonrap.sce new file mode 100755 index 000000000..63fc2abc8 --- /dev/null +++ b/50/DEPENDENCIES/newtonrap.sce @@ -0,0 +1,18 @@ + +function x=newton(x,f,fp) + R=5; + PE=10^-15; + maxval=10^4; + + for n=1:1:R + x=x-f(x)/fp(x); + + if abs(f(x))<=PE then break + end + if (abs(f(x))>maxval) then error('Solution diverges'); + abort + break + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/pivotgausselim.sci b/50/DEPENDENCIES/pivotgausselim.sci new file mode 100755 index 000000000..3fbd3402a --- /dev/null +++ b/50/DEPENDENCIES/pivotgausselim.sci @@ -0,0 +1,27 @@ +function [x]=pivotgausselim(A,b) + M=[A b]; + [ra,ca]=size(A); + [rb,cb]=size(b); + n=ra; + for p=1:1:n + for k=(p+1):1:n + if abs(M(k,p))>abs(M(p,p)) then + M({p,k},:)=M({k,p},:); + end + end + for i=p+1:1:n + m(i,p)=M(i,p)/M(p,p); + M(i,:)=M(i,:)-M(p,:)*m(i,p); + + end + end + a=M(1:n,1:n); + b=M(:,n+1); + for i = n:-1:1 + sumj=0 + for j=n:-1:i+1 + sumj = sumj + a(i,j)*x(j); + end; + x(i)=(b(i)-sumj)/a(i,i); + end +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/quadraticapprox.sci b/50/DEPENDENCIES/quadraticapprox.sci new file mode 100755 index 000000000..ad49267e9 --- /dev/null +++ b/50/DEPENDENCIES/quadraticapprox.sci @@ -0,0 +1,22 @@ +function [P]=quadraticapprox(x,f) + +n=length(x);m=length(f); +if m<>n then + error('linreg - Vectors x and f are not of the same length.'); + abort; +end; +s1=0; +s2=0; +for i=1:n + s1=s1+x(i)*f(i); + s2=s2+x(i)^2*f(i); +end +c0=det([sum(f) sum(x) sum(x^2);s1 sum(x^2) sum(x^3); s2 sum(x^3) sum(x^4)])/det([n sum(x) sum(x^2);sum(x) sum(x^2) sum(x^3); sum(x^2) sum(x^3) sum(x^4)]); + +c1=det([n sum(f) sum(x^2);sum(x) s1 sum(x^3); sum(x^2) s2 sum(x^4)])/det([n sum(x) sum(x^2);sum(x) sum(x^2) sum(x^3); sum(x^2) sum(x^3) sum(x^4)]); + +c2=det([n sum(x) sum(f);sum(x) sum(x^2) s1; sum(x^2) sum(x^3) s2])/det([n sum(x) sum(x^2);sum(x) sum(x^2) sum(x^3); sum(x^2) sum(x^3) sum(x^4)]); + +X=poly(0,"X"); +P=c2*X^2+c1*X+c0; +endfunction diff --git a/50/DEPENDENCIES/regulafalsi.sce b/50/DEPENDENCIES/regulafalsi.sce new file mode 100755 index 000000000..1b93787bf --- /dev/null +++ b/50/DEPENDENCIES/regulafalsi.sce @@ -0,0 +1,12 @@ +function [x]=regulafalsi(a,b,f) + N=100; + PE=10^-5; + for n=2:1:N + x=a-(a-b)*f(a)/(f(a)-f(b)); + if abs(f(x))<=PE then break; + elseif (f(a)*f(x)<0) then b=x; + else a=x; + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/regulafalsi4.sce b/50/DEPENDENCIES/regulafalsi4.sce new file mode 100755 index 000000000..5028fbdb0 --- /dev/null +++ b/50/DEPENDENCIES/regulafalsi4.sce @@ -0,0 +1,12 @@ +function [x]=regulafalsi4(a,b,f) + N=100; + PE=10^-5; + for n=2:1:N + x=a-(a-b)*f(a)/(f(a)-f(b)); + if abs(f(x))<=PE then break; + elseif (f(a)*f(x)<0) then b=x; + else a=x; + end + end + disp(n," no. of iterations =") +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/secant4.sce b/50/DEPENDENCIES/secant4.sce new file mode 100755 index 000000000..4cca19886 --- /dev/null +++ b/50/DEPENDENCIES/secant4.sce @@ -0,0 +1,12 @@ +function [x]=secant4(a,b,f) + N=4; // define max. no. iterations to be performed + PE=10^-4 // define tolerance for convergence + for n=1:1:N // initiating for loop + x=a-(a-b)*f(a)/(f(a)-f(b)); + if abs(f(x))<=PE then break; //checking for the required condition + else a=b; + b=x; + end + end + disp(n," no. of iterations =") // +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/secant64.sce b/50/DEPENDENCIES/secant64.sce new file mode 100755 index 000000000..b65a9f3d5 --- /dev/null +++ b/50/DEPENDENCIES/secant64.sce @@ -0,0 +1,12 @@ +function [x]=secant64(a,b,f,fp) + N=100; // define max. no. iterations to be performed + PE=10^-15 // define tolerance for convergence + for n=1:1:N // initiating for loop + x=(b*f(a)*fp(b)-a*f(b)*fp(a))/(f(a)*fp(b)-f(b)*fp(a)); + if abs(f(x))<=PE then break; //checking for the required condition + else a=b; + b=x; + end + end + disp(n," no. of iterations =") // +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/secant65.sce b/50/DEPENDENCIES/secant65.sce new file mode 100755 index 000000000..5d7d4e52c --- /dev/null +++ b/50/DEPENDENCIES/secant65.sce @@ -0,0 +1,13 @@ +function [x]=secant65(a,b,f) + deff('[y]=g(x)','y=-f(x)^2/(f(x-f(x))-f(x))'); + N=4; // define max. no. iterations to be performed + PE=10^-15 // define tolerance for convergence + for n=1:1:N // initiating for loop + x=a-(b-a)*g(a)/(g(b)-g(a)); + if abs(f(x))<=PE then break; //checking for the required condition + else a=b; + b=x; + end + end + disp(n," no. of iterations =") // +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/shooting.sci b/50/DEPENDENCIES/shooting.sci new file mode 100755 index 000000000..a9bf9bee6 --- /dev/null +++ b/50/DEPENDENCIES/shooting.sci @@ -0,0 +1,31 @@ +function [U] = shooting(ub,up,x,f) + +//Shooting method for a second order +//boundary value problem +//ub = [u0 u1] -> boundary conditions +//x = a vector showing the range of x +//f = function defining ODE, i.e., +// du/dx = f(x,u), u = [u(1);u(2)]. +//up = vector with range of du/dx at x=x0 +//xuTable = table for interpolating derivatives +//uderiv = derivative boundary condition + +n = length(up); +m = length(x); +y1 = zeros(up); + +for j = 1:n + u0 = [ub(1);up(j)]; + uu = ode(u0,x(1),x,f); + u1(j) = uu(1,m); +end; + +xuTable = [u1';up]; +uderiv = interpln(xuTable,ub(2)); +u0 = [ub(1);uderiv]; +u = ode(u0,x(1),x,f); +U=u'; + +endfunction + + diff --git a/50/DEPENDENCIES/simRK4.sci b/50/DEPENDENCIES/simRK4.sci new file mode 100755 index 000000000..d7dc4a913 --- /dev/null +++ b/50/DEPENDENCIES/simRK4.sci @@ -0,0 +1,27 @@ +function [u,v,t] = simRK4(u0,v0,t0,tn,h,f1,f2) + +// RK4 method solving simultanious ODE +// du/dt = f1(t,u,v), dv/dt = f2(t,u,v) with initial +//conditions u=u0,v=v0 at t=t0. The +//solution is obtained for t = [t0:h:tn] +//and returned in u,v + +umaxAllowed = 1e+100; + +t = [t0:h:tn]; u = zeros(t);v=zeros(t) ;n = length(u); u(1) = u0;v(1)=v0 + +for j = 1:n-1 + k11=h*f1(t(j),u(j),v(j)); + k21=h*f2(t(j),u(j),v(j)); + k12=h*f1(t(j)+h/2,u(j)+k11/2,v(j)+k21/2); + k22=h*f2(t(j)+h/2,u(j)+k11/2,v(j)+k21/2); + k13=h*f1(t(j)+h/2,u(j)+k12/2,v(j)+k22/2); + k23=h*f2(t(j)+h/2,u(j)+k12/2,v(j)+k22/2); + k14=h*f1(t(j)+h,u(j)+k13,v(j)+k23); + k24=h*f2(t(j)+h,u(j)+k13,v(j)+k23); + u(j+1) = u(j) + (1/6)*(k11+2*k12+2*k13+k14); + v(j+1) = v(j) + (1/6)*(k21+2*k22+2*k23+k24); + +end; + +endfunction diff --git a/50/DEPENDENCIES/sim_eulercauchy.sce b/50/DEPENDENCIES/sim_eulercauchy.sce new file mode 100755 index 000000000..b16af4509 --- /dev/null +++ b/50/DEPENDENCIES/sim_eulercauchy.sce @@ -0,0 +1,24 @@ +function [u,v,t] = simeulercauchy(u0,v0,t0,tn,h,f1,f2) + + +// du/dt = f1(t,u,v), dv/dt = f2(t,u,v) with initial +//conditions u=u0,v=v0 at t=t0. The +//solution is obtained for t = [t0:h:tn] +//and returned in u,v + + +umaxAllowed = 1e+100; + +t = [t0:h:tn]; u = zeros(t); n = length(u); u(1) = u0; + +for j = 1:n-1 + k11=h*f1(t(j),u(j),v(j)); + k21=h*f2(t(j),u(j),v(j)); + k12=h*f1(t(j)+h,u(j)+k11,v(j)+k21); + k22=h*f2(t(j)+h,u(j)+k11,v(j)+k21); + u(j+1) = u(j) + (k11+k12)/2; + v(j+1) = v(j) + (k21+k22)/2; + +end; + +endfunction diff --git a/50/DEPENDENCIES/simpson.sci b/50/DEPENDENCIES/simpson.sci new file mode 100755 index 000000000..da9992fd8 --- /dev/null +++ b/50/DEPENDENCIES/simpson.sci @@ -0,0 +1,3 @@ +function I=simpson(a,b,f) + I=((b-a)/6)*(f(a)+4*f((a+b)/2)+f(b)); +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/simpson13.sci b/50/DEPENDENCIES/simpson13.sci new file mode 100755 index 000000000..3e3c6b326 --- /dev/null +++ b/50/DEPENDENCIES/simpson13.sci @@ -0,0 +1,34 @@ +function [I] = simpson13(x,h,f) +//This function calculates the numerical integration of f(x)dx +//between limits x(1) and x(n) using Simpson's 1/3 rule +//Check that x and y have the same size (which must be an odd number) +//Also, the values of x must be equally spaced with spacing h +y=feval(x,f); +[nrx,ncx]=size(x) +[nrf,ncf]=size(y) +if ((nrx<>1)|(nrf<>1)) then + error('x or f, or both, not column vector(s)'); + abort; +end; +if ((ncx<>ncf)) then + error('x and f are not of the same length'); + abort; +end; +//check that the size of the lists xL and f is odd +if (modulo(ncx,2)==0) then + disp(ncx,"list size =") + error('list size must be an odd number'); + abort +end; +n = ncx; + +I = f(x(1)) + f(x(n)); +for j = 2:n-1 + if(modulo(j,2)==0) then + I = I + 4*f(x(j)); + else + I = I + 2*f(x(j)); + end; +end; +I = (h/3.0)*I +endfunction diff --git a/50/DEPENDENCIES/simpson38.sci b/50/DEPENDENCIES/simpson38.sci new file mode 100755 index 000000000..ca057dfa7 --- /dev/null +++ b/50/DEPENDENCIES/simpson38.sci @@ -0,0 +1,37 @@ +function [I] = simpson38(x,f) +//This function calculates the numerical integration of f(x)dx +//between limits x(1) and x(n) using Simpson's 3/8 rule +//Check that x and f have the same size (which must be of the form 3*i+1, +//where i is an integer number) +//Also, the values of x must be equally spaced with spacing h + +y=feval(x,f); +[nrx,ncx]=size(x) +[nrf,ncf]=size(y) +if ((nrx<>1)|(nrf<>1)) then + error('x or f, or both, not column vector(s)'); + abort; +end; +if ((ncx<>ncf)) then + error('x and f are not of the same length'); + abort; +end; +//check that the size of the lists xL and f is odd +if (modulo(ncx-1,3)<>0) then + disp(ncx,"list size =") + error('list size must be of the form 3*i+1, where i=integer'); + abort +end; +n = ncx; +xdiff = mtlb_diff(x); +h = xdiff(1,1); +I = f(x(1)) + f(x(n)); +for j = 2:n-1 + if(modulo(j-1,3)==0) then + I = I + 2*f(x(j)); + else + I = I + 3*f(x(j)); + end; +end; +I = (3.0/8.0)*h*I +endfunction diff --git a/50/DEPENDENCIES/straightlineapprox.sce b/50/DEPENDENCIES/straightlineapprox.sce new file mode 100755 index 000000000..9b8e50e19 --- /dev/null +++ b/50/DEPENDENCIES/straightlineapprox.sce @@ -0,0 +1,16 @@ +function [P]=straightlineapprox(x,f) + +n=length(x);m=length(f); +if m<>n then + error('linreg - Vectors x and f are not of the same length.'); + abort; +end; +s=0; +for i=1:n + s=s+x(i)*f(i); +end +c0=det([sum(f) sum(x); s sum(x^2)])/det([n sum(x);sum(x) sum(x^2)]); +c1=det([ n sum(f); sum(x) s])/det([n sum(x);sum(x) sum(x^2)]); +X=poly(0,"X"); +P=c1*X+c0; +endfunction diff --git a/50/DEPENDENCIES/taylor.sci b/50/DEPENDENCIES/taylor.sci new file mode 100755 index 000000000..ea5fdc776 --- /dev/null +++ b/50/DEPENDENCIES/taylor.sci @@ -0,0 +1,3 @@ +function u=taylor(t) + u=(t^1*U1)/fact(1)+(t^2*U2)/fact(2)+(t^3*U3)/fact(3)+(t^4*U4)/fact(4)+(t^5*U5)/fact(5)+(t^6*U6)/fact(6)+(t^7*U7)/fact(7)+(t^8*U8)/fact(8)+(t^9*U9)/fact(9)+(t^10*U10)/fact(10)+(t^11*U11)/fact(11) +endfunction \ No newline at end of file diff --git a/50/DEPENDENCIES/trapezoidal.sci b/50/DEPENDENCIES/trapezoidal.sci new file mode 100755 index 000000000..0ec3e7e49 --- /dev/null +++ b/50/DEPENDENCIES/trapezoidal.sci @@ -0,0 +1,7 @@ +// solves the definite integral by the trapezoidal rule, +// given the limits a,b and the function f, +// returns the integral value I + +function I =trapezoidal(a,b,f) + I=((b-a)/2)*(f(a)+f(b)); +endfunction \ No newline at end of file -- cgit