1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
|
/* {\hrulefill} *
{\bf Simpson}θ
{\ % beginning of TeX mode
\noindent
ʬ $\displaystyle\int_a^b f(x)\,dx$ ĤˡȤơ
ζɤȤ롣
$$\int_a^b f(x)\,dx\sim{h\over3}\,(y_0+4y_1+2y_2+4y_3+2y_4+4y_5+
\cdots +2y_{n-2}+4y_{n-1}+y_n)\ ,$$
\noindent
ǡ $n$ ϶Ȥ
$\displaystyle x_i=a+{i\over n}\,(b-a),\ y_i=f(x_i)\quad (i=0, 1, 2, \cdots, n)$
Ȥ롣
θξϡŪ䤵
äTaylor Ƴ
$$\int_{\xi-{1\over h}}^{\xi+{1\over h}}f(x)\,dx
\sim\int_{\xi-{1\over h}}^{\xi+{1\over h}}p(x)\,dx
={h\over3}\,\bigl(f(\xi-{1\over h})+4f(\xi)+f(\xi+{1\over h})\bigr)+O(h^5)$$
\noindent
$\,\displaystyle{n\over2}\,$ 碌뤳Ȥˤꡢ
θϾ롣
\special{epsfile=simpson.eps hscale=.7 vscale=.7 hoffset=25}
\vskip 13cm
\noindent
ǡ$p(x)$ ϼΤ褦ʾ룲¿༰ɽ:
$$
p\Bigl(\xi-{1\over h}\Bigr)=f\Bigl(\xi-{1\over h}\Bigr),
\ p(\xi)=f(\xi),
\ p\Bigl(\xi+{1\over h}\Bigr)=f\Bigl(\xi+{1\over h}\Bigr)\ .
$$
\noindent
{\tt A, B, F(X), N} Ѥơƿͤǿͼ¸ԤäƤߤ褦
% end of TeX mode }
* {\hrulefill} */
/* {\hrulefill\ simpson.c\ \hrulefill} */
#include <stdio.h>
#define A 0. /* ʬ $[a, b]$ {\hfill} */
#define B 1.
#define F(X) ((X)*(X)) /* ʬؿ $f(x)$ {\hfill} */
#define N 20 /* {\ $[a, b]$ ʬ $n$
\hfill} */
double simpson_rule() /* {\ Simpson θѤ
$\displaystyle\int_a^b f(x)\,dx$
ؿ \hfill} */
{
long i; /* long Ȥ {\hfill} */
double h, s, y[N + 1]; /* {\ å $h$, ʬ $s$,
ʬ $\displaystyle{i\over n}$
Ǥδؿ $y_i$ \hfill} */
h = 1. / (double) N; /* {\ $\displaystyle h={1\over n}$
\hfill} */
for (i = 0; i <= N; ++i) /* {\ $\displaystyle y_i
=f\bigl({i\over n}\bigr)
\quad (i=0, 1, 2, \cdots, n)$
\hfill} */
y[i] = F((double) i / (double) N);
/* ʲ {\hfill} */
/* $\displaystyle
s={h\over3}(y_0+4y_1+2y_2+4y_3+
\cdots +4y_{n-1}+y_n)$ {\hfill} */
/* η {\hfill} */
s = y[0];
for (i = 1; i < N; ++i)
{
if (i % 2 == 1) /* {\ {\tt i \% 2}= $i\ $
2 dzä; \hfill} */
s += 4. * y[i];
else
s += 2. * y[i];
}
s += y[N];
s *= (h / 3.);
return s; /* $s$ ֤ͤ {\hfill} */
}
main()
{
double s;
s = simpson_rule(); /* {\ ؿ
simpson\_rule() Ȥ \hfill} */
printf("%.16f\n", s); /* {\ $s$ ʲ16ɽ
\hfill} */
}
|
