Category: Numerical Analysis
Solution of the Laplace (Potential) Equation on a Two-Dimensional Square via Finite Differences
Solution of the One-Dimensional Heat Equation for a Rod Using Finite Differences by James Pate Williams, Jr. Created on Wednesday April 3, 2024
Undamped Mass-Spring Eigenvalue – Eigenvector Problem by James Pate Williams, Jr. (c) Monday April 1, 2024
We extend the results of the following website:
The five masses in the problem have a maximum value of 8. The six springs have a maximum value of 4 for their Hooke’s coefficients. The first 5 by 5 matrix is the inverse mass matrix, the second matrix is the Hooke’s coefficient 5 by 5 matrix, the third 5 by 5 matrix is the product of the inverse mass matrix times the Hooke’s coefficient matrix. The final row vector is the eigenvalue vector.
Lau Large Integer Calculator © March 27-28 by James Pate Williams, Jr. All Applicable Rights Reserved
Complex Number Calculator © March 25-26 by James Pate Williams, Jr.
New 1d Integration Results (c) March 24, 2024, by James Pate Williams, Jr.
I tested NUMAL’s integration function versus homegrown trapezoidal rule and Simpson’s rule. The second and third algorithms were closed (included both endpoints). There exist higher order Newton-Cotes integration formulas. I did not test the Gauss-Legendre and/or Gauss-Laguerre integration method(s). The trapezoidal rule can be improved if the derivative of the integrand is known and is easily calculated.
Triple Integration of Three Functions Using a Monte Carlo Algorithm and an Adaptive Quadrature Method (c) January 15 – 17, January 2024 by James Pate Williams, Jr.
Back in 2015 I translated a multiple integration FORTRAN subroutine to C# using switch statements to emulate conditional and unconditional go to statements. On January 15 – January 16, 2024, I translated my C# source code to vanilla C. Below is the FORTRAN subroutine’s website:
Here are some web pages with examples of triple integration:
3.4: Numerical Approximation of Multiple Integrals – Mathematics LibreTexts
3.3: Triple Integrals – Mathematics LibreTexts
https://tutorial.math.lamar.edu/Classes/CalcIII/IteratedIntegrals.aspx
https://tutorial.math.lamar.edu/Solutions/CalcIII/TripleIntegrals/Prob1.aspx
https://tutorial.math.lamar.edu/Problems/CalcIII/TripleIntegrals.aspx
I created a C test application to attempt to verify some of the results in online sources. For one example I also used my C# program.
== Menu ==
1 f(x, y, z) = 8 * x * y * z
2 f(x, y, z) = 4 * x * x * y - z * z * z
3 f(x, y, z) = x * y + z
5 Exit
Option: 1
f(x, y, z) = 8 * x * y * z
[0.000000, 1.000000] x [1.000000, 2.000000] x [2.000000, 3.000000]
N Integral -+Error % Error
10 0.754706 0.250631 94.968627
100 1.050263 0.126084 92.998246
1000 1.032125 0.036751 93.119165
10000 1.015068 0.011748 93.232880
100000 1.004080 0.003706 93.306131
1000000 1.000650 0.001169 93.329000
nQuadrature Integral Value and Percent Error: 15.000000 0.000000
== Menu ==
1 f(x, y, z) = 8 * x * y * z
2 f(x, y, z) = 4 * x * x * y - z * z * z
3 f(x, y, z) = x * y + z
5 Exit
Option:
== Menu ==
1 f(x, y, z) = 8 * x * y * z
2 f(x, y, z) = 4 * x * x * y - z * z * z
3 f(x, y, z) = x * y + z
5 Exit
Option: 2
f(x, y, z) = 4 * x * x * y - z * z * z
[2.000000, 3.000000] x [-1.000000, 4.000000] x [1.000000, 0.000000]
N Integral -+Error % Error
10 -12.087502 -3.929801 93.596026
100 -17.622371 -1.969516 90.663644
1000 -17.963837 -0.612528 90.482735
10000 -18.129559 -0.197880 90.394936
100000 -18.018231 -0.062690 90.453917
1000000 -17.939045 -0.019779 90.495870
nQuadrature Integral Value and Percent Error: 188.750000 200.000000
== Menu ==
1 f(x, y, z) = 8 * x * y * z
2 f(x, y, z) = 4 * x * x * y - z * z * z
3 f(x, y, z) = x * y + z
5 Exit
Option:
== Menu ==
1 f(x, y, z) = 8 * x * y * z
2 f(x, y, z) = 4 * x * x * y - z * z * z
3 f(x, y, z) = x * y + z
5 Exit
Option: 3
f(x, y, z) = x * y + z
[0.000000, 3.000000] x [0.000000, 2.000000] x [0.000000, 1.000000]
N Integral -+Error % Error
10 9.243858 2.121562 22.967849
100 12.062489 0.808150 0.520745
1000 12.096541 0.255335 0.804509
10000 12.094473 0.081030 0.787275
100000 12.051033 0.025704 0.425271
1000000 12.009588 0.008121 0.079900
nQuadrature Integral Value and Percent Error: 12.000000 0.000000
== Menu ==
1 f(x, y, z) = 8 * x * y * z
2 f(x, y, z) = 4 * x * x * y - z * z * z
3 f(x, y, z) = x * y + z
5 Exit
Option:
desired relative error 0.001
f# integral epsilon number err code
1 +1.4346639496E+000 1.295116E-003 257 0
2 +5.7531639665E-001 5.745470E-004 97 0
3 +2.1527578485E+000 1.793429E-003 45 0
4 +1.5998921741E+001 1.579802E-002 151 0
5 +1.8390615688E-001 6.510967E-005 97 0
6 -4.0003324629E+000 3.747971E-003 17 0
7 +8.6330831791E-001 8.563683E-004 45 0
8 +1.5000000000E+001 8.850520E-011 45 0
9 +1.8875000000E+002 1.109797E-009 45 0
10 +1.2000000000E+001 7.081136E-011 45 0
f# abs errror percent error
1 +9.7938837210E-005 +6.8261387483E-003
2 +4.7748257754E-005 +8.2987892410E-003
3 +6.1501589893E-004 +2.8576909005E-002
4 +1.0782593204E-003 +6.7391207526E-003
5 +9.9602348685E-007 +5.4159040098E-004
6 +3.3246288154E-004 +8.3115720386E-003
7 +2.6210055455E-004 +3.0369237392E-002
8 +8.8506979523E-011 +5.9004653015E-010
9 +1.1098109098E-009 +5.8797928998E-010
10 +7.0812689046E-011 +5.9010574205E-010
Numbers 8-10 in the preceding data correspond to the three three-dimensional functions in the menu illustrated above. The website claims option 2 in the Menu integral is -755 / 4 = -188.75. My calculation is the same magnitude but a positive sign.
Double Integration of Six Functions Using Three Different Quadrature Algorithms (c) January 15, 2024, by James Pate Williams, Jr. All Applicable Rights Reserved
I utilized six functions that have two independent variables. The integration methods were a Monte Carlo probabilistic procedure, TRICUB from “A Numerical Library in C for Scientists and Engineers” by H. T. Lau, PhD, and an algorithm from the ACM publication. I translated the ACM method from FORTRAN to C. See the following online references:
https://math.libretexts.org/Bookshelves/Calculus/Vector_Calculus_(Corral)/03%3A_Multiple_Integrals/3.04%3A_Numerical_Approximation_of_Multiple_Integrals
https://tutorial.math.lamar.edu/Classes/CalcIII/IteratedIntegrals.aspx
== Menu ==
1 8 * x + 6 * y
2 6 * x * y * y
3 2 * x - 4 * y * y * y
4 x * x * y * y + cos(pi * x) + sin(pi * y)
5 pow(2 * x + 3 * y, -2)
6 x * exp(x * y)
7 Exit
Option: 1
f(x, y) = 8x + 6y
[0.000000, 1.000000] x [0.000000, 2.000000]
N Integral +-Error % Error
10 18.398218 2.687822 8.008911
100 19.546981 0.838888 2.265096
1000 19.859442 0.266101 0.702789
10000 20.152819 0.082772 0.764093
100000 20.014598 0.026273 0.072988
1000000 20.013430 0.008323 0.067150
NUMAL Integral Value and Percent Error: 22.666667 13.333333
nQuadrature Integral Value and Percent Error: 19.999968 0.000160
== Menu ==
1 8 * x + 6 * y
2 6 * x * y * y
3 2 * x - 4 * y * y * y
4 x * x * y * y + cos(pi * x) + sin(pi * y)
5 pow(2 * x + 3 * y, -2)
6 x * exp(x * y)
7 Exit
Option:
== Menu ==
1 8 * x + 6 * y
2 6 * x * y * y
3 2 * x - 4 * y * y * y
4 x * x * y * y + cos(pi * x) + sin(pi * y)
5 pow(2 * x + 3 * y, -2)
6 x * exp(x * y)
7 Exit
Option: 2
f(x, y) = 6 * x * y * y
[1.000000, 2.000000] x [2.000000, 4.000000]
N Integral +-Error % Error
10 6.400786 2.133333 92.380017
100 7.528110 0.951714 91.037965
1000 7.951653 0.294436 90.533746
10000 8.083377 0.094959 90.376933
100000 8.001852 0.029859 90.473986
1000000 8.010457 0.009474 90.463742
NUMAL Integral Value and Percent Error: 204.800000 143.809524
nQuadrature Integral Value and Percent Error: 167.999730 99.999679
== Menu ==
1 8 * x + 6 * y
2 6 * x * y * y
3 2 * x - 4 * y * y * y
4 x * x * y * y + cos(pi * x) + sin(pi * y)
5 pow(2 * x + 3 * y, -2)
6 x * exp(x * y)
7 Exit
Option:
== Menu ==
1 8 * x + 6 * y
2 6 * x * y * y
3 2 * x - 4 * y * y * y
4 x * x * y * y + cos(pi * x) + sin(pi * y)
5 pow(2 * x + 3 * y, -2)
6 x * exp(x * y)
7 Exit
Option: 3
f(x, y) = 2 * x - 4 * y * y * y
[-5.000000, 4.000000] x [0.000000, 3.000000]
N Integral +-Error % Error
10 -645.214459 239.379239 14.654172
100 -436.170436 84.138767 42.305498
1000 -471.478924 25.861533 37.635063
10000 -486.066536 8.338088 35.705485
100000 -484.301805 2.644031 35.938915
1000000 -486.527596 0.839338 35.644498
NUMAL Integral Value and Percent Error: -765.000000 1.190476
nQuadrature Integral Value and Percent Error: -755.998734 0.000167
== Menu ==
1 8 * x + 6 * y
2 6 * x * y * y
3 2 * x - 4 * y * y * y
4 x * x * y * y + cos(pi * x) + sin(pi * y)
5 pow(2 * x + 3 * y, -2)
6 x * exp(x * y)
7 Exit
Option:
== Menu ==
1 8 * x + 6 * y
2 6 * x * y * y
3 2 * x - 4 * y * y * y
4 x * x * y * y + cos(pi * x) + sin(pi * y)
5 pow(2 * x + 3 * y, -2)
6 x * exp(x * y)
7 Exit
Option: 4
f(x, y) = x * x * y * y + cos(pi * x) + sin(pi * y)
[-2.000000, -1.000000] x [0.000000, 1.000000]
N Integral +-Error % Error
10 1.103184 0.173977 22.003269
100 0.708005 0.063869 49.942995
1000 0.757961 0.021546 46.411062
10000 0.744889 0.006866 47.335224
100000 0.748463 0.002171 47.082576
1000000 0.746892 0.000686 47.193609
NUMAL Integral Value and Percent Error: 1.449491 2.481177
nQuadrature Integral Value and Percent Error: 1.414395 0.000160
== Menu ==
1 8 * x + 6 * y
2 6 * x * y * y
3 2 * x - 4 * y * y * y
4 x * x * y * y + cos(pi * x) + sin(pi * y)
5 pow(2 * x + 3 * y, -2)
6 x * exp(x * y)
7 Exit
Option:
== Menu ==
1 8 * x + 6 * y
2 6 * x * y * y
3 2 * x - 4 * y * y * y
4 x * x * y * y + cos(pi * x) + sin(pi * y)
5 pow(2 * x + 3 * y, -2)
6 x * exp(x * y)
7 Exit
Option: 5
f(x, y) = pow(2 * x + 3 * y, -2)
[0.000000, 1.000000] x [1.000000, 2.000000]
N Integral +-Error % Error
10 0.555876 0.235032 1394.669315
100 0.752953 0.323309 1924.579264
1000 0.826769 0.275474 2123.058629
10000 0.927229 0.177631 2393.180482
100000 0.995193 0.097358 2575.927399
1000000 1.291253 0.127325 3371.987619
NUMAL Integral Value and Percent Error: 0.034125 8.243058
nQuadrature Integral Value and Percent Error: 0.037190 0.000639
== Menu ==
1 8 * x + 6 * y
2 6 * x * y * y
3 2 * x - 4 * y * y * y
4 x * x * y * y + cos(pi * x) + sin(pi * y)
5 pow(2 * x + 3 * y, -2)
6 x * exp(x * y)
7 Exit
Option:
== Menu ==
1 8 * x + 6 * y
2 6 * x * y * y
3 2 * x - 4 * y * y * y
4 x * x * y * y + cos(pi * x) + sin(pi * y)
5 pow(2 * x + 3 * y, -2)
6 x * exp(x * y)
7 Exit
Option: 6
f(x, y) = x * exp(x * y)
[0.000000, 1.000000] x [-1.000000, 2.000000]
N Integral +-Error % Error
10 3.217934 1.437426 inf
100 5.562365 0.873166 inf
1000 5.325378 0.225123 inf
10000 5.497069 0.075770 inf
100000 5.362919 0.023344 inf
1000000 5.371802 0.007410 inf
NUMAL Integral Value and Percent Error: 6.162907 inf
nQuadrature Integral Value and Percent Error: 2.562339 inf
Unknown option, please try again.== Menu ==
1 8 * x + 6 * y
2 6 * x * y * y
3 2 * x - 4 * y * y * y
4 x * x * y * y + cos(pi * x) + sin(pi * y)
5 pow(2 * x + 3 * y, -2)
6 x * exp(x * y)
7 Exit
Option:
Revisiting One-dimensional Definite Integration (c) January 15, 2024, by James Pate Williams, Jr.
These tests involve two integration methods from “A Numerical Library in C for Scientists and Engineers” by H. T. Lau, PhD, and two home-grown integration algorithms: Trapezoidal Rule and Simpson’s Rule. The second of the Lau algorithms performs the best.
Integral of f(x) = 10 / (x * x) for various a and b
INTEGRAL delivers:
-5.000000e+00 0 -5.000000e+00 -2 2
-7.499999e+00 0 -7.499999e+00 -4 1
-9.499999e+00 0 -9.499999e+00 -20 0
-9.999990e+00 1 -9.999990e+00 0 0
Integral of f(x) = sin(x) from 0 to pi
Delivers: 2.000000e+00 0
Approximation of the integral
x * x * exp(x) * sin(x) from 0 to pi
53.566456 0.000000
Steps Trapezoidal Simpson Err Trap Err Simp
16 52.835728 53.554295 1.364174 0.022724
32 53.383198 53.565685 0.342135 0.001460
48 53.484917 53.566315 0.152242 0.000285
64 53.520687 53.566395 0.085464 0.000135
80 53.537106 53.566418 0.054814 0.000093
96 53.546146 53.566395 0.037936 0.000135
Exact integral 53.566467
Numerical Explorations 