第七章 静电场与积分计算

如 |D| < , 则停止, 否则 i+=1, 继续。

稍加分析，我们会发现， Si+1 与 Si 的计算是 有联系的。 算法中 Si+1 的计算: 设计算 Si 时的步长为 H=(b-a)/2i, 则计算 Si+1 的步长为 H/2 . Si+1 的计算式 (n=2i+1+1):
1 3
( n 1 ) / 2 (n3) / 2


7.2变步长辛卜生积分
各种算法的积分运算，要想得到与实际积分足 够接近的数值结果，都必须取足够小的步长， 也即足够多的区间分割。 变步长辛卜生积分就是“逐次加倍等分区间的 个数，直至误差符合精度要求”的算法： 第1次, 把积分区间 2 等分, n=21+1=3, 结果记 S1 第2次, 把积分区间 4 等分, n=22+1=5, 结果记 S2 … 第i次, 把积分区间 2i 等分, n=2i+1, 结果记 Si 第i+1次, 把积分区间 2i+1 等分, n=2i+1+1, 结果 Si+1 ...
i 1
( H 是第 i 步时的步长)
2
R Pi [ f ( x1 ) f ( x n ) 2 RCi

1
f ( x1 2 k H )]
k 1

k 1
2
i 1
f ( x1 2 k H H )
所以有:
R Pi 1 R Pi 2 R C i
所以 Si+1 的计算过程为: 1. RP RP 2RC
或
V 1 3
( n 1) / 2 (n3) / 2
[ f ( x1 ) 4

i 1
f ( x2i ) 2

j 1
f ( x 2 j 1 ) f ( x n )] x
要求n为奇数, 即有(n-1)偶数个小条形.
梯形法和抛物线法都可以写成统一格式:
Vi Ci f ( xi )h
while (abs(D) > ee)|(H>H0) S1 = S2; RC1 = RC2; RP1 = RP2; RP2 = RP1+2*RC1; RC2 = sum(sin([H/2:H:(2^i-0.5)*H])); S2 = (RP2+4*RC2)*H/6; if abs(S2) < 1 D = S1-S2; else D = (S1-S2)/S2; end i = i+1; H = H/2; end v = S2;
2
2
数值积分原理: 把区间等分,求所有小条形面积的和,当x趋于0 时就是函数f(x)在整个区间上的积分值.
1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 -1
f(x)
x
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
矩形公式: b V f ( x)dx 对于积分: a 把区间[a,b]等分成n-1份,等分点记为x1, x2, …, xn 则等分间隔为: ba
x
矩形公式就是用矩形面积近似小条形的面积
n 1
xi
xi+1
每个小矩形面积: 矩形法计算公式:
S i f ( xi ) x
V
i 1, 2 , ..., n 1

n 1
f ( xi ) x
i 1
例: 求sin(x) 在[0, ]的积分.
function v = jx1(n) % 矩形法求sin(x)在[0,pi]上的积分 dx = pi/(n-1); x = 0:dx:pi; v = 0; for i = 1:n-1 v = v + sin(x(i))*dx; end
i 1 i
i
2. 3.
R C i 1

k 1
2
i
f ( x1 k H
H 2
)
S i 1
H 6
[ R Pi 1 4 R C i 1 ]
最后，为了防止假收敛,须引入H0, 当H< H0 时再 做收敛性判断.
function [v,i] = xps(H0,ee) % 变步长辛卜生法求sin(x)在[0,pi]上的积分 i = 1; RP1 = sin(0)+sin(pi); RC1 = sin(pi/2); H = pi/2; S1 = (RP1+4*RC1)*H/3; S2 = S1; RC2 = RC1; RP2 = RP1; D = 1;
k 1 k 1
)]
令:
R Pi 1 [ f ( x1 ) Fra Baidu bibliotekf ( x n ) 2 f ( x 1 k H )]
k 1
2 1
i
R C i 1

k 1
2
i
f ( x1 k H
H 2
)
则: 由于:
S i 1
H 6
[ R Pi 1 4 R C i 1 ],
S i 1
[ f ( x1 ) 4

k 1
f ( x2k ) 2

k 1
f ( x 2 k 1 ) f ( x n )]
H 2
其中:
x 2 k 1 x1 kH
x 2 k x1 k H
i
H 2
H 2
则:
S i 1 H 6
2 1 2
i
[ f ( x1 ) f ( x n ) 2 f ( x 1 k H ) 4 f ( x1 k H
这样得到的序列 {S1, S2, S3,…, Si, Si+1,…} 肯定是收敛于实际积分值的。然而具体 计算到哪一步结束，要看所需结果的精 度而定。 预先设置误差精度ε, 令

S i S i 1 , | S i 1 | 1 D S i S i 1 (相 对 误 差 ) , | S i 1 | 1 S i 1
抛物线形公式(n=2):
4 1 1 Vi f ( x i 1 ) f ( xi ) f ( x i 1 ) 2 x 6 6 6
例: 求sin(x) 在[0, ]的积分. function v = pwx(n) % 抛物线形法求sin(x)在[0,pi]上的积分 if mod(n,2)~=1 % 判断n是否是奇数. mod()求余数 函数. error('n should be odd'); end dx = pi/(n-1); x = 0:dx:pi; y = sin(x); v = (y(1)+y(n)+4*sum(y(2:2:n-1))+… 2*sum(y(3:2:n-2)))*dx/3;
抛物线形公式(辛卜生): 抛物线公式就是用抛物线代替f(x)进行积分
xi-1
xi
xi+1
梯形公式相当于用两点(xi-1,f(xi-1)), (xi,f(xi))上的 线形插值函数代替f(x)进行积分. 抛物线公式就是用三点(xi-1,f(xi-1)), (xi,f(xi)), (xi+1,f(xi+1)) 二次插值函数代替f(x)进行积分.
从而有：
Si

x i 1 x i 1
y ( x)dx
1 3
f ( x i 1 ) 4 f ( x i )
f ( x i 1 ) x
那么:
V 1 3 [ f ( x1 ) 4 f ( x 2 ) 2 f ( x 3 ) 4 f ( x 4 ) , ..., f ( x n )] x
function v = jx2(n) % 矩形法求sin(x)在[0,pi]上的积分 dx = pi/(n-1); x = 0:dx:pi; y = sin(x); v = sum(y(1:n-1))*dx; >> jx1(100) ans = 1.9998 >> jx2(100) ans =
1.8138
1 2 f ( x1 ) 1 2 f ( xn )
1 2
f ( x n )] x
V [

n 1
f ( x i ) ] x
i2
function v = tx(n) % 梯形法求sin(x)在[0,pi]上的积分 dx = pi/(n-1); x = 0:dx:pi; y = sin(x); v = (0.5*(y(1)+y(n))+sum(y(2:n-1)))*dx;
第七章 静电场与积分计算

静电场
y
( x0 , y0 )
线电荷电场
－A
A
x
V ( x0 , y0 )
1 4

A A
(x)
( x x0 ) y0
2 2
dx
取:
(x)
4
1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2
1, A 1则其等势线为:
-1.5
x xi 1 th, h ba 2 x
于是

b
a
y ( x)dx [
0
2
(t 1)(t 2) 2
f ( xi 1 ) t (t 2) f ( xi )
t (t 1) 2
f ( xi 1 )]hdt

h 3
[ f ( xi 1 ) 4 f ( xi ) f ( xi 1 )]
( x x i 1 )( x x i ) ( x i 1 x i 1 )( x i 1 x i )
f ( x i 1 )
即:
y(x) ( x x i )( x x i 1 ) 2(x)
2
f ( x i 1 )
( x x i 1 )( x x i 1 ) (x)
梯形公式: 梯形公式就是用梯形面积近似小条形的面积
xi
xi+1
Si
V
每个小梯形面积:
1 2
[ f ( x i ) f ( x i 1 )] x
1 2 [ f ( x i ) f ( x i 1 )] x
梯形法计算公式:

n 1
i 1
或:
即:
V [
1 2
f ( x1 ) f ( x 2 ) , ..., f ( x n 1 )
n i 0 n
Ci
n
1 (1)
n i
(t n i!(n i )!
0 j 0 j i
n
n
j )dt
称为牛顿-柯特斯公式, Cin 称为牛顿-柯特斯系数,与 f(x) 无关为常数. h 为区间范围. 梯形公式(n=1):
Vi [ 1 2 f ( xi ) 1 2 f ( x i 1 )] x
2
f ( xi )
( x x i 1 )( x x i ) 2(x)
2
f ( x i 1 )
取x
i 1
a, xi
ab 2
, xi 1 b ，在[a,b]上用
y ( x) f ( x) ，得

b
a
f ( x)dx y ( x)dx
a
b
为了便于计算，做变量代替
>> [v,i]=xps(0.1,0.01) v= 2.00000103336941 i= 5 >> [v,i]=xps(0.1,1e-10) v= 2.00000000000098 i= 10
等势线画法: function xdh % 等势线 V(x0,y0) = 1/4piee*int(l(x)/sqrt((xx0)^2+y^2),x,-1,1); % 线电荷，长为[-1,1],l(x) 为电荷线密度。 lpi=1/4piee*l(x). [X,Y] = meshgrid([-1.5:1/(4*pi):1.5]); n = 2001; dx = 2/(n-1); x = -1:dx:1; Z = zeros(size(X));
二次插值公式:
y(x) ( x x i )( x x i 1 ) ( x i 1 x i )( x i 1 x i 1 ) f ( x i 1 ) ( x x i 1 )( x x i 1 ) ( x i x i 1 )( x i x i 1 ) f ( xi )
-1
-0.5
0
0.5
1
1.5
取:
(x)
4
x ,A 1
2
1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2
-1.5
-1
-0.5
0
0.5
1
1.5

7.1 矩形梯形和抛物线形积分近似计算
V
定积分:

b a
f ( x ) d x F (b ) F ( a )
数值积分的原因: (1) f(x) 形式很复杂难于求其原函数F(x); 如: f ( x) ax bx c (2) f(x) 的原函数不能由初等函数表达; x 2 如: f ( x ) e , sin x (3) f(x) 没有给出表达式,只给出数据值.