北航數(shù)值分析第三次作業(yè)_圖文

1、數(shù)值分析第三次作業(yè)1. 算法的設(shè)計(jì)方案對(duì)Fredholm 積分方程，用迭代法進(jìn)行求解：u (x ' =A (u (x ，其中A (u (x =g (x -K (x , y u (y dy-11對(duì)于公式中的積分部分用數(shù)值積分方法。復(fù)化梯形積分法，取2601個(gè)節(jié)點(diǎn)，取迭代次數(shù)上限為50次。實(shí)際計(jì)算迭代次數(shù)為18次，最后算得誤差為r= 0.97E-10。復(fù)化Simpson 積分法，取迭代次數(shù)上限為50次，取2*41+1，即83個(gè)節(jié)點(diǎn)時(shí)能滿(mǎn)足精度要求。實(shí)際計(jì)算迭代次數(shù)為17次，最后的誤差為 r= 0.97E-10。Guass 積分法選擇的Gauss Legendre 法，取迭代次數(shù)上限為50次

2、，直接選擇8個(gè)節(jié)點(diǎn)，滿(mǎn)足精度要求。實(shí)際計(jì)算迭代次數(shù)為24次，最后算得誤差為r= 0.87E-10。2. 全部源程序 module integral implicit none contains!/復(fù)化梯形 subroutine trapezoid(m implicit none integer : i,j,k,mreal*8 : x(m+1,u(m+1 real*8 : sum,sum1,g,r,h real*8 : e=1.0e-10h=2./m do i=1,m+1x(i=-1.+(i-1*h end dou=0.02 do k=1,50 do i=1,m+1 sum1=0.g=dexp(

3、x(i*4.+(dexp(x(i+4.-dexp(-4.-x(i/(x(i+4. do j=2,m sum1=sum1+dexp(x(i*x(j*u(j end do sum=h/2.*(dexp(x(i*-1.*u(1+dexp(x(i*1.*u(m+1+2*sum1 u(i=g-sum end dor=h/2.*(dexp(x(1*4-u(1*2+(dexp(x(m+1*4-u(m+1*2 do i=2,mr=r+h*(dexp(x(i*4-u(i*2 end doif(dabs(r<=e exit end dowrite(*,* kopen(1,file="trapezo

4、id.txt" do i=1,m+1write(1,'(3(f18.12' x(i,u(i,dexp(x(i*4. end dowrite(1,'(4x,a2,e9.2' "r=",r close(1returnend subroutine trapezoid!/復(fù)化simpson subroutine simpson(m implicit none integer : i,j,k,mreal*8 : x(2*m+1,u(2*m+1 real*8 : sum,sum1,sum2,g,r,h real*8 : e=1.0e-10h=2

5、./(2.*m do i=1,2*m+1 x(i=-1.+(i-1*h end dou=0.02 do k=1,50do i=1,2*m+1 sum1=0. sum2=0.g=dexp(x(i*4.+(dexp(x(i+4.-dexp(-4.-x(i/(x(i+4. do j=1,msum1=sum1+dexp(x(i*x(2*j*u(2*j end do do j=1,m-1sum2=sum2+dexp(x(i*x(2*j+1*u(2*j+1sum=h/3.*(dexp(x(i*-1.*u(1+dexp(x(i*1.*u(2*m+1+4*sum1+2*sum2 u(i=g-sum end d

6、or=h/3.*(dexp(x(1*4-u(1*2+(dexp(x(2*m+1*4-u(2*m+1*2 do i=1,mr=r+4.*h/3.*(dexp(x(2*i*4-u(2*i*2 end do do i=1,m-1r=r+2.*h/3.*(dexp(x(2*i+1*4-u(2*i+1*2 end doif(dabs(r<=e exit end dowrite(*,* kopen(2,file="simpson.txt" do i=1,2*m+1write(2,'(3(f18.12' x(i,u(i,dexp(x(i*4. end dowrite

7、(2,'(4x,a2,e9.2' "r=",r close(2returnend subroutine simpson!/Gauss_Legendre法 subroutine Gauss implicit noneinteger,parameter : m=8 integer : i,j,kreal*8 : x(m,u(m,a(m real*8 : sum,g,r real*8 : e=1.0e-10data x /-0.9602898565,-0.7966664774,-0.5255324099,-0.1834346425,& 0.18343464

8、25,0.5255324099,0.7966664774,0.9602898565/data a /0.1012285363,0.2223810345,0.3137066459,0.3626837834,& 0.3626837834,0.3137066459,0.2223810345,0.1012285363/u=0.02 do k=1,50 do i=1,mg=dexp(x(i*4.+(dexp(x(i+4.-dexp(-4.-x(i/(x(i+4. do j=1,msum=sum+dexp(x(i*x(j*u(j*a(j end do u(i=g-sum end do r=0.do

9、 i=1,mr=r+a(i*(dexp(x(i*4-u(i*2 end doif(dabs(r<=e exit end dowrite(*,* kopen(3,file="Gauss.txt" do i=1,mwrite(3,'(3(f18.12' x(i,u(i,dexp(x(i*4. end dowrite(3,'(4x,a2,e9.2' "r=",r close(3returnend subroutine Gaussend module!/主程序 program main use integral implicit noneinteger : code1=2600 integer : code2=41call trapezoid(code1 call simpson(code2 call Gaussend program3. 各種積分方法的節(jié)點(diǎn)和數(shù)值解（由于數(shù)據(jù)太多，在打印時(shí)用了較計(jì)算時(shí)少的有效數(shù)字）復(fù)化梯形法 復(fù)化Simpson 法 4. 各方法所得曲線(xiàn) （由于所取節(jié)點(diǎn)太多，且精度高，所以圖中很難看出各曲線(xiàn)的區(qū)別。 ） 5. 不