длина дуги кривой в Фортране

Программа должна вычислить длину кривой ƒ=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

fortran fortran95
person Nikos Ellinakis    schedule 17.11.2015    source источник
comment
Недостаточно времени, чтобы просмотреть код, но я бы порекомендовал хотя бы определить индекс цикла i как целое число и вычислить только x, объявленный как double precision, из i.   -  person Vladimir F    schedule 17.11.2015
comment
ошибка порядка операций, (5.3*20/x*abs(x2-x1)) должно быть (5.3*20/(x*abs(x2-x1))) (конечно, определение функции или предварительное вычисление локального x позволяет избежать такой ошибки)   -  person agentp    schedule 17.11.2015

Ответы (1)


arrow_upward
2
arrow_downward

Всякий раз, когда вы сталкиваетесь с проблемами, связанными с функцией, пожалуйста, создайте функцию. Это определенно сделает жизнь бесконечно проще!

Это поможет вам в отладке, так как вы сможете сразу выяснить, действительно ли ваша функция выполняет правильные вычисления.

Есть много проблем с вашим кодом, вы смешиваете целые числа и вещественные числа, никогда не делайте этого. В конце концов, это создаст вам проблемы. В идеале все должно быть определено с помощью 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

Подводя итог, кодирование не всегда связано с прямой реализацией, но приведенный выше код ясен, и вы легко найдете любые ошибки, поскольку каждая строка имеет не более нескольких операций, пожалуйста, попробуйте придерживаться этой методологии, она будет обязательно поможет потом...

person zeroth    schedule 17.11.2015
comment
zeroth, Большое спасибо за вашу помощь, я очень ценю это, хотя у меня есть еще один вопрос. Будет ли программа работать нормально, если я смогу прочитать x1,x2 с клавиатуры? - person Nikos Ellinakis; 17.11.2015
comment
@NikosEllinakis Уверяю вас, учиться, пытаясь, — это самый способ научиться программированию. - person zeroth; 17.11.2015