算出的y2结果错了，请求帮忙！

马媛

算出来的y1是正确的，但是算y2时算错了。找不到原因了。请大家帮我看看！

fcode

请给出 ess.dat 文件

vvt

OCR 工具，嘿嘿

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

Program qmzy11
  Implicit None
  Integer , parameter :: N = 500
  Integer :: i
  Real(Kind=8) :: x0(N), y0(N) , x1 , y1 , x2 , y2
  Open (11, File='ess.dat')
  Open (12, File='dt.dat')
  Do i = 1, 50
    Read (11, *) x0(i), y0(i)
  End Do
  x1 = 0.2655D0
  Call lagr(x0, y0, N, x1, y1)
  x2 = 0.3886D0
  Call lagr(x0, y0, N, x2, y2)
  Write (12, 1) x1, x2
  Write (*, 1) y1, y2
  1 Format (7X, 2E15.7)
End Program qmzy11

Subroutine lagr(x0, y0, n, x, y)
  Implicit None
  Integer :: i , j , n
  Real(Kind=8) :: x0(n), y0(n) , p , x , y
  y = 0.D0
  Do i = 1, n
    p = 1.0D0
    Do j = 1, n
      If (i==j) cycle
      p = p*(x-x0(j))/(x0(i)-x0(j))
    End Do
    y = y + p*y0(i)
  End Do
End Subroutine lagr