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

第4章附录
4.2.2 复化求积分
例题4.2.5计算程序
//simp.c//
# include<stdio.h>
# include<math.h>
# 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))
write(2,'(2x,2f10.6)') x(i),f(i)
end do
a = 1.;
b = 4.;
c = 1.
do i = 1,n-1
r(i) = 3.0*(f(i+1)-f(i-1))/dx
end do
r(1) = r(1)-fp(0); r(n-1) = r(n-1)-fp(n)
call tridag(a,b,c,r,u,n-1) ! 注意三对角矩阵是n-1维u(0) = fp(0); u(n) = fp(n)
do i = 0,n
write(5,'(2x,2f10.6)') x(i),u(i)
end do
end
subroutine tridag(a,b,c,r,u,n)
parameter (nmax=500)
integer n,j
real(8) a(0:n),b(0:n),c(0:n),r(0:n),u(0:n)
real(8) bet,gam(nmax)
if(b(1).eq.0.)pause 'tridag: rewrite equations'
bet=b(1)
u(1)=r(1)/bet
do j=2,n
gam(j)=c(j-1)/bet
bet=b(j)-a(j)*gam(j)
if(bet.eq.0.)pause 'tridag failed'
u(j)=(r(j)-a(j)*u(j-1))/bet
end do
do j=n-1,1,-1
u(j)=u(j)-gam(j+1)*u(j+1)
end do
return
end
!!!Romberg.for
program romberg
external f
double precision f,a,b,s
a=0.0
b=1.0
eps=0.000001
call romb(a,b,f,eps,s,n)
write(*,10) s,n
10 format(1x,'s=',d15.8,1x,'n=',i6)
end program romberg
function f(x)
double precision f,x
! f=log(1+x)/(1+x*x)
! f=x/(4+x*x)
f=4./(1.+x*x)
return
end
subroutine romb(a,b,f,eps,t,n)
dimension y(10)
double precision a,b,f,t,y,h,p,s,q
h=b-a
y(1)=h*(f(a)+f(b))/2.0
m=1
n=1
10 p=0.0
do 20 i=0,n-1
20 p=p+f(a+(i+0.5)*h)
p=(y(1)+h*p)/2.0
s=1.0
do 30 k=1,m
s=4*s
q=(s*p-y(k))/(s-1)
y(k)=p
p=q
30 continue
if ((abs(q-y(m)).ge.eps).and.(m.le.9)) then
m=m+1
y(m)=q
n=n+n
h=h/2.0
goto 10
end if
t=q
return
end。