工程分析程序设计 上机作业（五）模块化编程上机目的：掌握内部例程、主程序、外部例程、模块等功能的使用方法。

1、 有一个六边形，求其面积。

为求面积，作了3条辅助线。

如图所示：（提示，三角形面积=()()()s s a s b s c ---，其中2a b c s ++=，a 、b 、c 为三个边长）。

要求用内部函数来计算每个三角形的面积。

l 1=10 l 2=30l 3=16 l 4=13l 5=21l 6=14l 7=20l 8=36 l9=28program mainimplicit nonereal :: A(9)=(/10,30,16,13,21,14,20,36,28/)integer ireal ss=area(A(3),A(4),A(7))+area(A(2),A(7),A(8))+& area(A(1),A(8),A(9))+area(A(5),A(6),A(9))print*,'The area is',scontainsfunction area(x,y,z)implicit nonereal x,y,zreal p,areap=(x+y+z)/2area=sqrt(p*(p-x)*(p-y)*(p-z))end functionend2、编写一个子例程子程序，SUM(S，T，N1，N2)。

把整型数N1到N2进行求和，并把求和的结果放置到S，把整型数N1到N2进行求积，并把求积的结果放置到T。

并用这个子程序来计算：y=++++++++++-(1234)(345678)(3*4*5*6)(1*2*3)program mainimplicit noneinteger A(8)integer i,b,c,d,e,f,g,m,n,yb=0;d=0;g=1;n=1do i=1,8A(i)=iend docall sum(b,c,1,4)call sum(d,e,3,8)call sum(f,g,3,6)call sum(m,n,1,3)y=b+d+g-nprint*,'y=(1+2+3+4)+(3+4+5+6+7+8)+(3*4*5*6)-(1*2*3)=',ycontainssubroutine sum(s,t,N1,N2)integer s,t,N1,N2,ido i=N1,N2s=s+A(i)t=t*A(i)end doend subroutineend3、编写函数子程序GDC求两个数的最大公约数。

求最大公约数的算法如下：把两个数中大的那个数作为被除数，两数相除得到一个余数。

把余数去除除数得到新一轮的余数。

不断重复这一过程直到余数为0，这时的除数就是两个数的最大公约数。

调用此函数，求1260，198，72三个数的最大公约数。

program mainimplicit noneinteger a,b,c,m,n,tempprint*,'请输入三个数：'print*,'a b c'read*,a,b,cif (a<b) thentemp=aa=bb=tempend ifcall GDC(a,b,m)if (c<m) thentemp=cc=mm=tempend ifcall GDC(m,c,n)print*,'最大公约数为',ncontainssubroutine GDC(x,y,z)integer x,y,z,kdo while(mod(x,y)/=0)k=mod(x,y)x=yy=kend doz=yend subroutineend4、 编写一个模块程序，提供以下服务：定义出常量π、e 。

定义出子程序，实现求和21n i n=∑、求阶乘n !。

并在主程序中计算如下结果：从键盘上输入整数n 、实型数A 、R 、R 0，求21!n i n n =⎛⎫ ⎪⎝⎭∑（实型）和0202n n R R An R e R R π⎛⎫- ⎪⎝⎭⎛⎫ ⎪⎝⎭（实型）。

module mytaskimplicit nonereal,parameter :: PI=3.1415926real,parameter :: e =2.7182818containsfunction sub(n)integer i,n,subsub=0do i=1,nsub=sub+i*iend doend function subfunction fact(n)integer i,nreal factfact=1do i=1,nfact=fact*iend doend function factend module mytaskprogram mainuse mytaskimplicit noneinteger nreal A,R,R0real x,yprint*,'请输入依次 n,A,R,R0 :'read*,n,A,R,R0x=fact(n)/sub(n)y=(A*n/(2*PI*R*R))*(R/R0)**n*e**(-(R/R0)**n)print*,'n!/∑n^2(i=1,n)=',xprint*,'(A*n/(2*PI*R*R))*(R/R0)**n*e**(-(R/R0)**n)',yend program5、 编写外部函数，计算...!7!5!3sin 753+-+-=x x x x x 的值，直到最后一项的绝对值小于10-6为止。

再编写主程序，从键盘读入x ，调用该外部函数并输出sinx 的计算结果。

注意不能用Fortran 的标准函数SIN(X)。

program mainimplicit nonereal x,pexternal subprint*,'请输入角度 x :'read*,xp=x/180*3.14call sub(p)end programsubroutine sub(k)implicit nonereal k,yinteger iy=0i=1do while (abs(k**i/(fact(i)))>1.0e-6)y=y+(-1)**((i-1)/2)*k**i/(fact(i))i=i+2end doprint*,'sin(x) = ',ycontainsfunction fact(m)integer m,j,factfact=1do j=1,mfact=fact*jend doend function factend subroutine sub6、 编写程序，用欧拉法数值求解一阶微分方程。

数值解法的基本思想：用差分方程代替微分方程，然后在若干个离散点上逐点求解差分方程，得到各离散点x 0, x 1, x 2, … 等处函数(),()dy f x y x dx=的近似值y 0, y 1, y 2, … 。

其中各离散点x 0, x 1, x 2, …之间的距离称为步长h 。

欧拉法的基本思想是，一阶微分dy dx 可用向前差分代替，即11'()n n y y y h+≈-。

带入微分方程，可得 1(,)n n n n y y f x y h+-≈ 因此，1*(,)n n n n y y h f x y +≈+，其中n n x x h -=+1。

当给定x 0, y 0和步长h 后，即可按下面步骤求得数值解：0000(,)x f x y y ⎫⇒⇒⎬⎭1111,(,)y f x y x ⎫⇒⎬⎭()2222,n y f x y y x ⎫⋅⋅⋅⇒⎬⎭请按上述方法编写程序用Euler 法求解微分方程22dy y x dx=-，当x 0=0时，y 0=1.0。

取h =0.1，试求出当x i =0.1,0.2,0.3,0.4….,1.0时的y i 值。

program mainimplicit nonereal x(0:10)real y(0:10)real p(0:10)integer(2) i,jy(0)=1.0; x(0)=0.0p(0)=f(x(0),y(0))do i=1,10x(i)=x(i-1)+0.1end dodo j=1,10y(j)=y(j-1)+0.1*p(j-1)p(j)=f(x(j),y(j))print*,'x(',j,')=',x(j),'y(',j,')=',y(j)end docontainsfunction f(a,b)real a,b,ff=b*b-a*aend functionend。