Try porting the "FORTRAN77 Numerical Computing Programming" program to C and Python (Part 2)
Previous article Porting the "FORTRAN77 Numerical Programming" program to C and Python (Part 1)
From the code of FORTRAN77 Numerical Calculation Programming, 1.2.
Cumulative error in trapezoidal law
Cumulative rounding error
Quoted from the end of P.5 to the front of the P.6 program
Integral as an example of cumulative rounding error
I = \int_{0}^{1}\frac{1}{1+x^2}dx = \frac{\pi}{4} = 0.7853982\tag{1.13}
>, ** Trapezoidal rule ** while gradually reducing the step width $ h $
>```math
I_h = h\biggr[\frac{1}{2}f(a)+\sum_{j=1}^{n-1}f(a+jh)+\frac{1}{2}f(b)\biggr], \quad h=\frac{b-a}{n}\tag{1.14}
a=0, \quad b=1, \quad f(x)=\frac{1}{1+x^2} \tag{1.15}
Try to integrate numerically with>.
Fortran code
A program that performs numerical integration of (1.14) and (1.15).
strap.f
* Accumulation of round off error
PROGRAM STRAP
PARAMETER (ONE = 1.0)
*
EV = ATAN(ONE)
*
WRITE (*,2000)
2000 FORMAT (' ---- TRAP ----')
*
DO 10 K = 1, 13
*
N = 2**K
H = 1.0 / N
*
S = 0.5 * (1.0 + 0.5)
DO 20 I = 1, N - 1
S = S + 1 / (1 + (H * I)**2)
20 CONTINUE
S = H * S
*
ERR = S - EV
IF (ERR .NE. 0.0) THEN
ALERR = ALOG10 (ABS(ERR))
ELSE
ALERR = -9.9
END IF
*
WRITE (*,2001) N, H, S, ERR, ALERR
2001 FORMAT (' N=',I6,' H=',1PE9.2,' S=',E13.6,
$ ' ERR=',E8.1,' L(ERR)=',0PF4.1)
*
10 CONTINUE
*
END
The execution result is
---- TRAP ----
N= 2 H= 5.00E-01 S= 7.750000E-01 ERR=-1.0E-02 L(ERR)=-2.0
N= 4 H= 2.50E-01 S= 7.827941E-01 ERR=-2.6E-03 L(ERR)=-2.6
N= 8 H= 1.25E-01 S= 7.847471E-01 ERR=-6.5E-04 L(ERR)=-3.2
N= 16 H= 6.25E-02 S= 7.852354E-01 ERR=-1.6E-04 L(ERR)=-3.8
N= 32 H= 3.12E-02 S= 7.853574E-01 ERR=-4.1E-05 L(ERR)=-4.4
N= 64 H= 1.56E-02 S= 7.853879E-01 ERR=-1.0E-05 L(ERR)=-5.0
N= 128 H= 7.81E-03 S= 7.853957E-01 ERR=-2.4E-06 L(ERR)=-5.6
N= 256 H= 3.91E-03 S= 7.853976E-01 ERR=-5.4E-07 L(ERR)=-6.3
N= 512 H= 1.95E-03 S= 7.853978E-01 ERR=-4.2E-07 L(ERR)=-6.4
N= 1024 H= 9.77E-04 S= 7.853981E-01 ERR=-6.0E-08 L(ERR)=-7.2
N= 2048 H= 4.88E-04 S= 7.853982E-01 ERR= 6.0E-08 L(ERR)=-7.2
N= 4096 H= 2.44E-04 S= 7.853983E-01 ERR= 1.2E-07 L(ERR)=-6.9
N= 8192 H= 1.22E-04 S= 7.853981E-01 ERR=-6.0E-08 L(ERR)=-7.2
Will be. Oh, the output of ERR and L (ERR) is different from the book from around $ N = 128 $ ...
Understanding the code
Why is the name of this program STRAP in the first place? Aside from the question, what was the integral? Let's take a second look and check what we are writing in the book and the code. Numerical integration of the integral formula written in the book by the trapezoidal rule is to put a lot of trapezoids in the area of the graph represented by the original integration and combine the sum of the areas of the trapezoids into the integration graph. To make it an approximation of the area of.
So, what the program is doing is executing the formulas (1.14) and (1.15) while making $ h $ smaller and smaller. The theme is how much it differs from the answer required in (1.13).
I drew a flow chart.
First, what we are doing in * 1 is $ arctan 1.0 $, which uses the inverse trigonometric function arctangent. Expressing the arc tangent as a definite integral,
arctan\,x = \int_0^x \frac 1 {z^2 + 1}\,dz
It becomes the formula. If $ x $ is $ 1 $ and $ z $ is $ x $, it will be the same as the formula in (1.13), so the same value as (1.13), that is, $ I $ will be assigned to EV here. ..
Next, in * 2, $ 2 ^ K $ is substituted for $ N $, and $ N $ is the formula (1.14) for
And for * 3, this was originally expressed as
Then, the first expression in (1.14)
I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}f(a)+\sum_{j=1}^{n-1}f(a+j\frac{1}{n})+\frac{1}{2}f(b)\biggr]
And if you replace (1.15) with $ a = 0 $ and $ b = 1 $
I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}f(0)+\sum_{j=1}^{n-1}f(0+\frac{j}{n})+\frac{1}{2}f(1)\biggr]
And because
I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}\times\frac{1}{1+0^2}+\sum_{j=1}^{n-1}\frac{1}{1+\big(\frac{j}{n}\big)^2} +\frac{1}{2}\times\frac{1}{1+1^2}\biggr]
From
I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}+\sum_{j=1}^{n-1}\frac{n}{1+\frac{j^2}{n^2}}+\frac{1}{4}\biggr]\tag{r1}
further
I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}+\sum_{j=1}^{n-1}\frac{n^2}{n^2+j^2}+\frac{1}{4}\biggr]
so,
I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{3}{4}+\sum_{j=1}^{n-1}\frac{n^2}{n^2+j^2}\biggr]\tag{r2}
Can be organized.
- 4 substitutes $ 0.5 \ times (1.0 + 0.5) $ for $ S
, which is substituted for $ \ frac {1} {2} \ times \ big (\ frac {2} {2} + \ frac {1} {2} \ big)$ Think of it as $ \ frac {1} {2} \ times \ frac {3} {2} = \ frac {3} {4} $$ And it is the value that appears in the formula (r2).
The counter I in loop 2 is $$ j $ written in
And in * 5,
So, loop 2 is (1.14)
\frac{1}{2}f(a)+\sum_{j=1}^{n-1}f(a+jh)+\frac{1}{2}f(b)
It will be the one that is calculating the part of.
In * 6, the first equation of (1.14) is completed, and in * 7, the difference from the equation of (1.13), that is, the difference between the calculation result originally obtained by integral calculation and the approximate value using the trapezoidal law. I'm calculating.
In * 8, when the difference obtained in * 7 is not equal to 0, the absolute value of the difference is calculated (ʻABS (ERR) ), and the common logarithm of that value is calculated (ʻALOG10 ()). ..
The common logarithm is the value of 
This graph shows that the cumulative error decreases in proportion to $ h ^ 2 $, but it does not decrease as the error approaches the machine epsilon.
Difference from books
By comparison, the value of $ S $ has increased slightly from $ N = 128 $. The point where the decrease in cumulative error stops is a little later. The computer used in the book (dare not written as a computer) is MV4000 AOS / VS, so it's a 16-bit OS, isn't it? So, it's similar to the graph I put out, the graph calculated using the double precision floating point number that is given as a reference in the book. It's a little different just because it's similar, but in the book environment, floating-point arithmetic is 6-digit rounding arithmetic. So, I guess this is due to the result of truncation operation on a 64-bit OS (the environment is 32-bit). Please let me know if you are an expert. Postscript: I really had an expert (@cure_honey) tell me in the comments. Thank you very much.
C code
strap.c
#include <stdio.h>
#include <math.h>
int main(void)
{
const float ONE = 1.0f;
float ev = atanf(ONE);
int k, i;
int n;
float h, s, err, alerr;
printf("%14.7E\n", ev);
printf("--- TRAP ---\n");
for(k=1; k<=13; k++)
{
n = pow(2,k);
h = 1.0f / n;
s = 0.5f * (1.0f + 0.5f);
for(i=1; i<=n-1; i++)
{
s = s + 1 / (1+ powf(h*i,2));
}
s = h * s;
err = s - ev;
if(err != 0.0f)
{
alerr = log10(fabsf(err));
}
else
{
alerr = -9.9f;
}
printf("N=%6d H=%9.2E S=%13.6E ERR=%8.1E L(ERR)=%4.1F\n",n, h, s, err, alerr);
}
return 0;
}
The execution result does not change.
Python code
strap.py
import numpy as np
import matplotlib.pyplot as plt
ONE = np.array((1.0),dtype=np.float32)
ev = np.arctan(ONE,dtype=np.float32)
h, s, err, alerr = np.array([0.0, 0.0, 0.0, 0.0],dtype=np.float32)
#For calculation 1.0, 0.5, 0.0, -9.9 single precision floating point number
tmp = np.array([1.0, 0.5, 0.0, -9.9],dtype=np.float32)
#List for graph drawing
x = []
y = []
print("--- TRAP ---")
for k in range(13):
n = np.power(2,k+1)
h = tmp[0] / n.astype(np.float32)
s = tmp[1] * (tmp[0] + tmp[1])
for i in range(n-1):
s = s + tmp[0] / (tmp[0] + np.square(h*(i+1),dtype=np.float32))
s = s * h
err = s - ev
if err != tmp[2]:
alerr = np.log10(np.abs(err))
else:
alerr = tmp[3]
print("N=%6d H=%9.2E S=%13.6E ERR=%8.1E L(ERR)=%4.1F" % (n, h, s, err, alerr))
#Variable set for graph
x.append(k+1)
y.append(alerr)
#Graph drawing
fig, ax = plt.subplots()
ax.scatter(x,y)
ax.set_xlabel("n (1/2^n)")
ax.set_ylabel("l (10^l)")
plt.show()
The output result is the same. Also, I put out a graph.
Summary
I think the Python code is going to be cleaner. And I think it's a good idea to scatter numpy for single precision. Do you want to double-precision the Fortran code ... No, but then the execution result cannot be compared with the book. Muu.
Recommended Posts