计算程序

c name：qun-yue.for，无界承压含水层干扰井群法预测开采降深场，有潜水越流进入；当有界时，则把虚井同实井一样代入即可（但其Q却有正负之分）

c 计算公式：第一越流系统计算公式——汉图士-雅可布公式

c 即：s（j）＝∑Qi*F（ui，r/B）/（4πT），u（i）＝S*r（i，j）**2/（4Tt）

c B＝sqrt（T*m’/k’）r/B＝r/sqrt（T*m’/k’）越流系数＝k’/m’

c 其中越流井函数F（u，r/B），本程序自动查表求其值；

c 均质，各向同性s：降深（m），S：储水系数；m：预测点个数n：抽水井个数（包括实虚）ihv：预测时段个数；T：承压水全区T值（m*m/d），W：承压水全区S值，bk1m：弱透水层越流系数，即K’/m’（1/d）

parameter（m＝4，n＝2，ihv＝40，ihv1＝34）

double precision u1

c 抽水井的x、y坐标（单位：m）预测点的x、y坐标（单位：m），预测点可以是网格状剖分节点加非节点的预测点

dimension xn（n），yn（n），r（n，m）

dimension xm（m），ym（m）

dimension q（n），s（m，ihv），tt（ihv）

dimension bu（40），brb（30），bf（40，30）

character file3（m）*8，fileb*8，file2（m）*4

cmm：越流井函数F（u，r/B）表列数，nn：越流井函数F（u，r/B）表行数

mm＝28

nn＝38

open（1，file＝’biaof.txt’）

read（1，*）（brb（j），j＝1，mm）

do 210 i＝1，nn

read（1，*）bu（i），（bf（i，j），j＝1，mm）

210 continue

close（1）

c 抽水井点（n眼）坐标

open（1，file＝’xy-jing.txt’）

do 9 i＝1，n

read（1，*）xn（i），yn（i）

9 continue

close（1）

c 预测点（m个）坐标

open（1，file＝’xy-yuce.txt’）

do 10 i＝1，m

read（1，*）xm（i），ym（i）

10 continue

close（1）

open（1，file＝’kmwb.txt’）

read（1，*）ak，am，w，bk，bm

close（1）

T＝ak*am

bk1m＝bk/bm

b＝sqrt（T/bk1m）

open（1，file＝’q.txt’余烂）

read（1，*）（q（i），i＝1，n）

close（1）

c’yucehao.txt’存放m个字符串，每串4字符，为预测点输出预测的历时水位曲线标注字符用

open（1，file＝’yucehao.txt’）

read（1，42）（file2（i），i＝1，m）

close（1）

42 format（5a4）

c’yucewen.txt’存放m个字符串，为＂预测号.dat＂，每串8字符，为预测点输出预测的历时水位曲线准备文件名

open（1，file＝’yucewen.txt’）

read（1，44）（file3（i），i＝1，m）

close（1）

44 format（5a8）

c’tl1’-预测掘帆时间系列：累计时间（d）

open（1，file＝’tl1.txt’）

read（1，*）（tt（i），i＝1，ihv）

close（1）

tt0＝0

c 各预测点j至抽水井点i之距r（m）

do 30 j＝1，m

s（j，0）＝0

do 20 i＝1，n

r（i，j）＝sqrt（（xn（i）-xm（竖散漏j））**2+（yn（i）-ym（j））**2）

if（r（i，j）.eq.0）r（i，j）＝0.15

20 continue

30 continue

close（1）

open（1，file＝’sids.txt’）

do 60 k＝1，ihv1

write（*，*）’k＝’，k，’t＝’，tt（k）

write（1，70）k，tt（k）

70 format（/2x，’K＝’，i3，’t＝’，f15.6，’日’）

c 下一大段对m个预测点逐点计算1～n眼抽水井在k时刻对每个预测点水位降之和s（j，k）

do 50 j＝1，m

s（j，k）＝0

c 计算r/B，及k时刻i～j泰斯井函数自变量u1

do 40 i＝1，n

u1＝r（i，j）**2*w/（4.0*T*tt（k））

rb＝r（i，j）/b

c 限制u1超出8

if（u1.gt.8）u1＝8

c 限制r/B超出9

if（rb.gt.9）rb＝9

c 越流井函数F（u，r/b），程序自查表求之

c 程序自查表，先查列，查计算值rb位于F表的哪两列brb（jj-1）～brb（jj）之间

do 220 jj＝2，mm

jx＝jj

if（rb.ge.brb（jj-1）.and.rb.le.brb（jj））goto 230

220 continue

c 计算出F表的brb（jj-1）～brb（jj）两列之间的列距dx

230 dx＝brb（jx）-brb（jx-1）

c 计算值rb到F表brb（jj-1）列之距dx1、到F表brb（jj）列之距dx2

dx1＝rb-brb（jx-1）

dx2＝brb（jx）-rb

c 程序自查表，再查行，查计算值u1位于F表的哪两行bu（ii-1）～bu（ii）之间

do 240 ii＝2，nn

iy＝ii

if（u1.ge.bu（ii-1）.and.u1.le.bu（ii））goto 250

240 continue

c 计算出F表的bu（ii-1）～bu（ii）两行之间的行距dy

250 dy＝bu（iy）-bu（iy-1）

c 计算值u1到F表bu（ii-1）行之距dy1、到F表bu（ii）行之距dy2

dy1＝u1-bu（iy-1）

dy2＝bu（iy）-u1

c 找出计算值（u1，rb）在F表中最近四点的F值（即bf），其中，f11为左上角，f12右上角，f21左下角，f22右下角

f11＝bf（iy-1，jx-1）

f12＝bf（iy-1，jx）

f21＝bf（iy，jx-1）

f22＝bf（iy，jx）

f1＝f11+dx1/dx*（f12-f11）

f2＝f21+dx1/dx*（f22-f21）

f＝f1+dy1/dy*（f2-f1）

c 用越流泰斯公式计算k时刻第i眼抽水井对第j个预测点的水位降ss

ss＝q（i）/（4.0*3.1415926*T）*f

c 把第i眼（i＝1～n）抽水井在k时刻对第j预测点水位降ss累加到j点的迭加降深s（j，k）中

s（j，k）＝s（j，k）+ss

40 continue

50 continue

do 100 j＝1，m 100 write（1，80）xm（j），ym（j），-s（j，k），j

60 continue

close（1）

80 format（2f10.1，f9.3，i4）

do 300 j＝1，m

fileb＝file3（j）

open（1，file＝fileb）

write（1，310）tt0，s（j，0），file2（j）

do 320 k＝1，ihv1

write（1，380）tt（k），-s（j，k）

320 continue

300 close（1）

380 format（f10.1，f8.2）

310 format（f10.1，f8.2，1x，a4）

cq＝0

do 90 i＝1，n

90 cq＝cq+q（i）

open（1，file＝’qq.txt’）

write（1，140）（i，q（i），i＝1，n）

write（1，*）

write（1，*）’总采量＝’，cq，’（m3/d）’

close（1）

140 format（5（i4，f8.0））

open（1，file＝’okk.dat’）

do 290 j＝1，m

290 write（1，80）xm（j），ym（j），-s（j，ihv1），j

close（1）

stop

end

用二分法求方程x*x-x-1=0的正根，要求精确到小数点后四位。（matlab）

l1 计算公式

f(ak)*f(bk)<0

bk-ak=1/2k-1*(b-a)

a1<=a2<=… <=ak<=…,b1<=b2<=…<=bk<=…。

l2 算法分析

设f(x)∈C[a,b],且f(a)f(b)<0，在[a,b]间寻找实根，记[a1,b1]=[a,b]，取x1=(a1+b1)/2,若f(x1)=0,则x1是f(x)=0的根，f(x1)f(a1)>0,则a2=x1,b2=b1,否则a2=a1,b2=x1。得到[a2,b2]满足：f(a2)f(b2)<0,b2-a2=(b1-a1)/2=(b-a)/2,a2>=a1,b2<=b1。以[a2,b2]取代[a1,b1]，继续以上过程，直到精度达到要求为止。

l3 源程序

function f1=fun(x)

f1=x-cos(x)

function [x,k]=erfen(a,b,s) %a,b为根区间，s为精度

a=0b=1s=1e-4k=0

while abs(a-b)>s

x=(a+b)/2

if fun(a)*fun(x)<0

b=x

else

a=x

end

k=k+1

end

x=(a+b)/2 %x为方程的解

k % k为计算次数

实验结果讨论和分析

本题使用二告纯拍分法得到的x=0.7391，满足基本要求，题目要求精确裤悄到小数点袜羡后四位，告诉了本题二分法得应达到得精确度；计算次数为14，二分法收敛性很好，收敛速度不快。

---