Fortran中曲线的弧长

Question

程序必须计算ƒ=3.1*x^2-5.3/x在x=1/2和x=3/2之间的曲线长度。该长度应计算为n线段的总和，以n=1开始并以n=20结束.

我真的找不到为什么我得到的结果是错误的。例如，如果 x1=1/2 和 x2=3/2 得到 110，而我应该得到 13 我给你下面的代码：

program pr2_ex2
  implicit none

  integer::x
  double precision::dy,dx !dy=the height of the linear part & dx=the lenght of the linear part
  double precision::x1,x2,s !f=f(x) the function,x=the values which can be given to f
  double precision::length 

  print*,"Please enter the function's starting point"
  read*,x1

  print*,"Please enter the function's ending  point"
  read*,x2

  length = 0
  s = 0

  do x = 2, 21
    dx = ((x*abs(x2-x1)-(x-1)*abs(x2-x1))/(20))
    dy = (3.1*(x*abs(x2-x1)/20)**2-(5.3*20/x*abs(x2-x1)))-(3.1*((x-1)*abs(x2-x1)/20)**2-(5.3*20/(x-1)*abs(x2-x1)))

    length = sqrt((dx**2)+(dy**2))
    s = length+s
   end do

   print*,s

end program

Answer 1

每当您遇到涉及函数的问题时，请做创建一个函数。它一定会让生活变得更加轻松！

它将帮助您进行调试，因为您可以立即确定您的函数是否确实进行了正确的计算。

你的代码有很多问题，你把整数和实数混合在一起，永远不要那样做。它最终会给你带来问题。理想情况下，应使用 selected_real_kind 和相关的精度名称 0._dp 或 0._sp 或任何您将命名的精度来定义所有内容。

这是一个计算段长度的代码，结果刚好在 13 以上。

program pr2_ex2
  implicit none
  integer, parameter :: dp = selected_real_kind(p=15)

  integer :: i
  double precision :: dy, dx, dxsq
  double precision :: x1, x2, s, x
  double precision :: length

  integer :: N 

  ! perhaps read in N as well?
  N = 20

  print*,"Please enter the function's starting point"
  !  read*,x1
  x1 = 0.5_dp
  print*,"Please enter the function's ending  point"
  ! read*,x2
  x2 = 1.5_dp

  ! are you allowed to have x2 < x1, if so abort?
  ! dx is always the same, so why calculate it all the time?
  dx = abs(x2 - x1) / real(N,dp)
  ! we need not calculate this all the time
  dxsq = dx ** 2

  ! Total length
  length = 0._dp
  do i = 1 , N 

     ! starting x of this segment
     x = x1 + dx * (i-1)

     ! Get current delta-y
     dy = f(x + dx) - f(x)

     ! calculate segment vector length
     s = sqrt(dxsq + dy ** 2)

     ! sum total length
     length = length + s

  end do

  print*,length

contains

  function f(x) result(y)
    double precision, intent(in) :: x
    double precision :: y

    y = 3.1_dp * x ** 2 - 5.3_dp / x

  end function f

end program pr2_ex2

归根结底，编码并不总是直接实现，但上面的代码很清楚，你会很容易发现引入的任何错误，因为每一行只有几个操作，请尝试并坚持这种方法，以后肯定会对你有所帮助...