第4章 数值微分与积分(附录)

第4章附录

4.2.2 复化求积分

例题4.2.5计算程序

//simp.c//

# include

# include

# define f(x) 4./(1+(x)*(x))

void main(void)

{ float a = 0., b = 1., s, h;

int n = 100, i;

h = (b-a)/n;

s = f(a)+f(b);

for(i=1;i<=n-1;i++)

{ if(i%2==0)

s = s+2.0*f(a+i*h);

else

s=s+4.*f(a+i*h);

}

s = s*h/3;

printf("%10.5f\n",s);

}

====================================

4.2.3 变步长求积分公式和龙贝格求积分公式

例题4.2.6计算程序

!!!!Trapezia.for!!!

program trapezia

external f

real(8) f,a,b,s

a=0.0; b=1.0; eps=1.e-6

call trap(a,b,f,eps,s,n)

write(*,10) s,n

10 format(1x,'s=',d15.6,3x,'n=',i5)

end

function f(x)

real(8) x

f=exp(-x*x)

end

subroutine trap(a,b,f,eps,t,n)

real(8) a,b,f,t,fa,fb,h,t0,s,x

fa=f(a); fb=f(b)

n=1; h=b-a

t0=h*(fa+fb)/2.0

5 s=0.0

do 10 k=0,n-1

x=a+(k+0.5)*h

s=s+f(x)

10 continue

t=(t0+h*s)/2.0

if (abs(t-t0).ge.eps) then

t0=t

n=n+n

h=h/2.0

goto 5

end if

return

end

%%% demo_aTrapInt.m %%%

function demo_aTrapInt

clc;clear; format long;

[T nsub] = aTrapInt(@f01,0,1,0.000001)

end

function [T nsub]= aTrapInt(f,a,b,eps)

tol = 1; nsub = 1;

inall = 0;

T = 0.5*(b-a)*(f(a)+f(b));

while tol > eps

T0 = T;

nsub = 2*nsub;

n = nsub + 1; % total number of nodes h = (b-a)/nsub; % stepsize

x = a:h:b; % divide the interval inall = inall + sum(f(x(2:2:n-1)));

T = 0.5*h * (f(a)+2*inall+f(b));

tol = abs(T-T0);

end

end

%%% demo_aSimpInt.m%%%

function demo_aSimpInt

clc;clear; format long;

[S nsub] = aSimpInt(@f01,0,1,0.000001)

end

function [S nsub] = aSimpInt(f,a,b,eps)

nsub = 2;

even = f((a+b)/2);

odd = 0;

S = (b-a)*(f(a) + 4*even + 2*odd + f(b))/6;

inall = even + odd;

tol = 1;

while tol > eps

S0 = S;

nsub = 2*nsub;

n = nsub + 1; % total number of points

h = (b-a)/nsub; % stepsize

x = a:h:b; % divide the interval

even = sum(f(x(2:2:n-1))); % 偶节点

odd = inall; % 间隔取半前的全部内节点inall 在间隔取半时全部变为奇节点S = (h/3)*( f(a) + 4*even +2*odd + f(b));

inall = even + odd;

tol = abs(S - S0);

end

====================================

例题4.2.7计算程序

!!!Simpson.f90

program simpson ! f=sin(x),fp=cos(x)

parameter(pi=3.1415926,n=64)

real(8) a(0:n),b(0:n),c(0:n),r(0:n),u(0:n)

real(8) x(0:n),f(0:n),fp(0:n),dx,d

integer i

open(2,file='exa4_1_3_old.txt')

open(5,file='exa4_1_3.txt')

dx = 2.0*pi/n

d = 3.0/dx

fp(0) = 1.0 ! fp(0) = cos(0) =1

fp(n) = 1.0 ! fp(n) = cos(2pi)=1

do i = 0,n

x(i) = i*dx; f(i) = sin(x(i))