Nikos Ellinakis
Nikos Ellinakis

Reputation: 1

arc length of a curve in Fortran

The program must calculate the length of the curve of ƒ=3.1*x^2-5.3/x between x=1/2 and x=3/2.The length should be calculated as the sum of n line segments starting with n=1 and ending with n=20.

I really can't find why the result I'm getting is wrong. For example if x1=1/2 and x2=3/2 then I get 110 when I should be getting 13 I am giving you the code below:

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

Upvotes: 0

Views: 383

Answers (1)

nickpapior
nickpapior

Reputation: 782

Whenever you face problems involving a function, then please, do create a function. It will definitely make life endlessly more easy!

It will help you debug as you can immediately figure out whether your function actually does the correct calculation.

There are many issues with your code, you are mixing integers and real, never do that. It will give you problems eventually. Ideally everything should be defined using selected_real_kind and the associated precision designation 0._dp or 0._sp or whatever you will name your precisions.

Here is a code which calculates the segment length, and yields just above 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

Bottom line, coding is not always about being direct in implementation, but the above code is clear, and you will easily find any bugs introduced as each line has no more than a few operations, please do try and adhere to this methodology, it will most definitely help you later...

Upvotes: 2

Related Questions