计算方法的作业——请网友点评一下 - Fortran Coder

[Fortran] 纯文本查看 复制代码

!/*************************************************************************

!   > File Name: kt0301_2.f90

!   > Author: 青羽100P 

!   > Mail: 3443988515@qq.com 

!   > Created Time: 2017年10月11日 星期三 18时03分27秒

! ************************************************************************/



      program kt0301_2

!拉格朗日主程序代码

         implicit none

         integer::pa(17),pa_1(4),i,z,error

         real::tem(17),tem_1(4),chads

         !使用数组存放数值

         open(10,file='sample_unformatted_sequential.dat',f&

             &orm='unformatted',access='sequential',status='old')

         i=1

         do

           read(10,iostat=error),tem(i),pa(i)

           if(error/=0)exit

           i=i+1

         end do

         close(10)

         pa(1)=1000

         pa_1=(/pa(1),pa(4),pa(7),pa(10)/)

         tem_1=(/tem(1),tem(4),tem(7),tem(10)/)

         !挑出1000，700，400，200处的数值作为已知点

         !-------------------------

         open(11,file="chazhishuju.txt")

         do z=1,17

             call lar(4,pa_1,tem_1,pa(z),chads)

             write(11,"(i2,f10.4,i5,f10.4,f10.4)"),z,tem(z),pa(z),chads,chads-tem(z)

         end do

      end program

      include 'lar.f90'

[Fortran] 纯文本查看 复制代码

!/************************************************************************

 !  > File Name: lar.f90

 !  > Author: 青羽100P 

 !  > Mail: 3443988515@qq.com 

 !  > Created Time: 2017年10月10日 星期二 21时59分40秒

 ! *************拉格朗日子程序代码***********************************************************/



      subroutine lar(n,x,y,cha,ds)

         implicit none

          integer::n,x(n),cha,i,j

          !n表示数组长度，x横坐标，cha表示要差值x的位置

          real::y(n),larlc,ds

          ds=0

          larlc=1

          do j=1,n

          larlc=1

          !保证不受上一次运算结果影响

              do i=1,n

                 if(i==j)then

                     larlc=larlc*1.0

                 else

                     larlc=larlc*((cha-x(i)*1.0)/(x(j)*1.0-x(i)*1.0))

                 end if

                 !此处if 块，实现连乘算式，当算到自己是等于1就可以了，也可以直接用cycle

                     !每个算式后 *1.0 保证实型运算，否则会被省略，结果全为零 

              end do

              ds=ds+larlc*y(j)

          end do

          print*,ds

      end subroutine