Solving Equation by Iteration

[Andrew_F_1]

Hello,

I wrote the following code to solve an equation by iteration but I never get a reasonable answer.

do f=0.0000,180.0000
y= (Thita_ow*pi/180d0) + (f*pi/180d0) - pi - ((cos(Thita_ow*pi/180d0)*cos((Thita_ow-Half_Angle)*pi/180d0))/(sin(Half_Angle*pi/180))) &
   + ((2*cos(Thita_ow*pi/180d0)-cos(f*pi/180d0))*((cos((f+Half_Angle)*pi/180d0))/(sin(Half_Angle*pi/180d0))))
if (y==0.0) exit
end do

! Creating the TxT file
open(unit=3 , file='test.txt')
write(3,*) f, y
close(unit=3)

Any help ?

[mecej4]

Your description leaves out many details, but a guess can be made as to the failure of the program. The probability is zero that the loop will be exited because y becomes exactly zero, unless the values of the known variables (Thita_ow, Half_Angle, etc.) have serendipitous values that make y = 0.

When the loop is exited, f will have the value 181, and y will be the value calculated with f=180.0. Describe the problem mathematically and show the complete program, and you may receive more helpful suggestions.

[John_Campbell]

You may consider the changes I have made below to reformat your computation (probably thwarted by IDZ).
As Mecej4 indicates, the test of y==0.0 is a bad test and you need to consider the accuracy you wish to achieve. I have provided a value y_tol for this purpose.
Alternatively, given that your loop is looking for the nearest degree of F, you might alternatively find the angle F with the minimum absolute value of y.
Running the sample shows that your test of y==0.0 does not work and the testing of each value f is a bad approach. A better solution, assuming that the function is continuous is to use a binary search, where given two values of f, f1=0 and f2 = 180, then test f = (f1+f2)/2 and if the sign of y(f) = sign of y(f1) then replace f1 with f, else f2 with f. Keep going until abs(f2-f1) < f_tol. I'll leave that to you.

My suggestion for an initial change is

!  test of loop
!
  integer*4 f
  real*8    Thita_ow, Half_Angle, y, pi
  real*8    Thita_ow_rad, Half_Angle_rad, y_tol, f_rad
!
! Variables Provided
  Thita_ow    = 75
  Half_Angle  = 90
!
  pi             = 4 * atan (1.0d0)
  y_tol          = 4 * epsilon (1.0d0)
!
  Thita_ow_rad   = Thita_ow   * pi/180d0
  Half_Angle_rad = Half_Angle * pi/180d0
!
  do f=0,180
    f_rad = f * pi/180d0
!
!    y = (Thita_ow*pi/180d0)       &
!      + (f*pi/180d0)              &
!      - pi                        &
!      - ( cos(Thita_ow*pi/180d0) * cos((Thita_ow-Half_Angle)*pi/180d0) / sin(Half_Angle*pi/180) ) &
!      + ( ( 2*cos(Thita_ow*pi/180d0)-cos(f*pi/180d0) ) * cos((f+Half_Angle)*pi/180d0) / sin(Half_Angle*pi/180d0) )
!
    y = Thita_ow_rad              &
      + f_rad                     &
      - pi                        &
      - ( cos(Thita_ow_rad)                * cos(Thita_ow_rad-Half_Angle_rad) / sin(Half_Angle_rad) )       &
      + ( (2*cos(Thita_ow_rad)-cos(f_rad)) * cos(f_rad+Half_Angle_rad)        / sin(Half_Angle_rad) )
!
    if (abs(y) < y_tol) exit
     write (*,*) f, y
  end do
!
! Creating the TxT file
  open (unit=13 , file='test.txt')
  write (13,*) f, f_rad, y
  close (unit=13)

end

John

[John_Campbell]

Sorry to harp on about this, but if I use cut and paste from a notepad, I really do not understand why IDZ does not cope with this and provide a reasonable representation of what was posted. Most windows text files have lines terminated with <cr> <lf> and to be able to test for this and reproduce this in the same way that most other display programs can do should not be a big ask. Surely, can't you test for <cr> <lf> and reproduce the text accordingly ?

To fail to fix this is not attending to the needs of we windows users of this forum and treating we clients with contempt.

After all, this is a windows user forum which provides over 50% of all posts !!

John

[John_Campbell]

The attached program solves for y small for f in the range 0 to 180 deg by using a binary search.
I have assumed values for Thita_ow and Half_Angle and that there is a single zero value between Y(0) and Y(180) and that the function Y is continuous. (Math 101). You could expand on the code for when this is not the case.

John

[andrew_4619]

@John > I have shared your pain on posting code but have had success recently , the tip is:

In the posting window there is a menu bar, use the icon with { } and 'code' below it. This pops a new window and in a pull down option you pick the type of code (ie fortran), and then past into the popup box. This seems to  fix most of my issues. before I was using the fortran tags which does not work well.

   real*8    Thita_ow, Half_Angle, f, f1, f2, y, y1, y2, f_tol, y_tol, y_rad
   integer*4 i
   external  y_rad
!
! Variables Provided
   Thita_ow   = 35
   Half_Angle = 90
   y_tol      = 4 * epsilon (1.0d0)
   f_tol      = y_tol
!
! Initial guess for F1, F2; y1 x y2 must be < 0
   f1 = 0
   f2 = 180
!
   y1 = y_rad ( f1, Thita_ow, Half_Angle)
   y2 = y_rad ( f2, Thita_ow, Half_Angle)
!
   if ( y1*y2 >= 0) then
      write (*,*) 'BAD initial choice for f1,f2', f1,f2, y1,y2
      stop
   end if
!  
   do i = 1,200
!
! Use a binary search
     f  = ( f1+f2 ) / 2
     y  = y_rad ( f, Thita_ow, Half_Angle)
!
     if ( y*y1 > 0) then
        f1 = f
        y1 = y
     else
        f2 = f
        y2 = y
     end if
!
     if ( abs(f2-f1) < f_tol ) then
        write (*,*) f,y, ' f1,f2 too close'
        exit
     else if (abs(y) < y_tol) then
        write (*,*) f,y, ' y within tollerance'
        exit
     end if
      write (*,*) i, f, y
   end do
 
! Creating the TxT file
   open (unit=13 , file='test.txt')
   write (13,*) f, y
   close (unit=13)
 
 end
 
 real*8 function y_rad ( f_deg, Thita_ow_deg, Half_Angle_deg )
!
   real*8    f_deg, Thita_ow_deg, Half_Angle_deg
   real*8    f_rad, Thita_ow_rad, Half_Angle_rad, pi
!
   pi             = 4 * atan (1.0d0)
!
   f_rad          = f_deg          * pi/180d0
   Thita_ow_rad   = Thita_ow_deg   * pi/180d0
   Half_Angle_rad = Half_Angle_deg * pi/180d0
!
   y_rad = Thita_ow_rad              &
         + f_rad                     &
         - pi                        &
         - ( cos(Thita_ow_rad)                * cos(Thita_ow_rad-Half_Angle_rad) / sin(Half_Angle_rad) )       &
         + ( (2*cos(Thita_ow_rad)-cos(f_rad)) * cos(f_rad+Half_Angle_rad)        / sin(Half_Angle_rad) )
!
 end function y_rad

[John_Campbell]

@app4619,

Thanks very much for the advice. I tried it and it appears to work well.

I only hope others become more aware of the {...} code option.

John

[jimdempseyatthecove]

Another gripe issue is you used to be able to end a line with Ctrl-Enter
And have the next line abutted up tight against the prior line.
During composing of the Comment it appears so.
However, after submitting it does not.
Sometimes I have lists of things to enter the I do not wish to use bullets.
This post made using Ctrl-Enter

Jim Dempsey

[jimdempseyatthecove]

Yes, you can use the {...} and choose Plain Text

But then
it appears in a goofy box like this

Out of line - out of context.

Jim Dempsey

[Andrew_F_1]

Guys,

I modified my code and it is working now.

do f=0.0,180.0,0.00001
y= (Thita_ow*pi/180d0) + (f*pi/180d0) - pi - ((cos(Thita_ow*pi/180d0)*cos((Thita_ow-Half_Angle)*pi/180d0))/(sin(Half_Angle*pi/180))) &
   + ((2*cos(Thita_ow*pi/180d0)-cos(f*pi/180d0))*((cos((f+Half_Angle)*pi/180d0))/(sin(Half_Angle*pi/180d0))))
if (y>0.0000001) exit
end do

Knowing that:

Half_Angle= 30

Thita_ow= 180

The answers I get are f=   66.79241   and  y=  2.8448312E-07

This code satisfy my needs but it is kinda memory consuming. Thanks all, appreciate it!

[mecej4]

Your algorithm, although it works for now, is neither robust nor efficient -- it may be compared to searching for the lost MH370 airplane using a microscope; if the expression whose zero is sought is a decreasing function and the first trial value yields a positive value of y, the loop will terminate prematurely and give you an incorrect "answer".

A good method such as Mueller's method can produce the result after evaluating the function less than ten times, whereas your simple search algorithm evaluates the function several million times!

[John_Campbell]

I looked at a change to the binary search as provided in Quote #6 and replaced the mid point by a linear extrapolation rule.

 f  = ( f1*y2 - f2*y1 ) / (y2 - y1)        ! linear search

By searching in the range -180:180 with half_angle=30 and thita_ow = 180; by retaining the sign strategy for replacing f1 or f2, this converged in 17 iterations, but if instead of the sign test, I replace the larger of y1 or y2 then it does not converge. This compares with 50 iterations with binary search for the same precision and less than 10 for Mueller's method.

There is something to be said for a robust strategy, which the binary search has, as the sign replacement rule ensures there is always a zero value between f1 and f2.

I would suspect that Mueller's method of quadratic extrapolation would work well for well behaved functions, but it does not take too much for a function to misbehave, and that is not including the problem where the function becomes discontinuous.

John

[Kevin_D_Intel]

Jim, John - I passed along your Forum usage feedback to the Forum Development/Support team.

[jimdempseyatthecove]

John,

A thought came to me just now, I haven't fully fleshed it out, (ignore it if you wish).

The methods listed above are scalar in nature, and has a sizable computation per iteration incorporating sin and cos.

What you might consider doing is rework the code to use vectors by way of the IPP library ippCos_... and ippSin_... functions.

On a CPU with SSEn.nn you could produce two results per iteration, AVX 4 results, MIC 8 results.

The binary search could then become 3-way, 5-way, or 9-way search, thus halving, quartering, eighth-ing the search time with the additional overhead of migrating the best choice to the other lanes of the vector.

On the SSE you choose he 1/3 and 2/3 division points as opposed to the 1/2 point of the formal binary search. Perform the expression on the two halves of the vector at the same time. Then a little scalar code (or swizzle code) to make the selection of the winner to populate the other lane(s) of the vector, then make your termination test on y and break the loop when converges.

AVX you choose the 1/5, 2/5, 3/5 and 4/5 points. 

MIC/AVX512 1/9, 2/9, ...

Jim Dempsey

[DavidWhite]

Further to John's comments in Quote 12, instead of binary searching, Fibonnaci search is more efficient.  Given the tolerance required on the final solutoin, the number of iterations required is known before starting.  There is an overhead, of calculating the Fibonacci numbers before the search begins, but it is a robust and efficicent method, which does not rely on the function behaving itself.

Regards,

David

[Anthony_Richards]

It appears that the variable is F as you are trying values of F from 0 to 180 and that F represents an angle in degrees.

Since you are restricting your F range, you are assuming that there is a solution for F within this range which gives Y=0. This is not necessarily going to be the case, given that the other quantities Theta and Halfangle introduce another two double infinities of possible values.

Therefore you need an algorithm that will allow for there (a) not being a solution in the range 0->180 and (b) not being any solution at all.

A standard way of approaching such a non-linear equation as you have is to use a Taylor series expansion to second order in small quantities about a starting value for F and solve for a small change in F that will move F towards a value which allows Y to approach zero.
Thus if F0 is a first guess then you solve 0=Y(F0) +dF0*(dY/dF) +0.5*(dF0)^2*(d2Y/dF2)  for dF0 where the coefficients are the first and second differential coefficients of your function taken with respect to the variable F (and evaluated at F=F0). This is a simple quadratic equation A*dF0^2+B*dF0+C=0 with standard solutions which may or may not be real, depending on the values for the coefficients.

Your function Y is a differentiable function of the variable F and so you can form functions for computing dY/dF and d2Y/dF2. It is then a case of iterative (recursive) solution: Start with F0, compute dF0, compute F1=F0+dF0, compute dF1, compute F2=F1+dF1 and so on, while testing the absolute magnitude of successive Y(Fn,Theta, Halfangle) values against a small number near zero until it falls below that value, when you then have found a fairly accurate solution (accuracy depending on the magnitude of the small number). At each stage your function for solving the quadratic must return a flag value that indicates when no solution is possible.

I have written a simple version of your equation which converges to 10^-10 after four iterations from a starting guess of F=90. I simplified your equation to Y(f,theta,phi)=sin(phi)*(theta+F-pi) -(cos(theta)*cos(theta-phi)) + (2*cos(theta)-cos(F))*cos(F+phi) which I think is a correct reading where I have used phi instead of your Halfangle. I set Theta=75 degrees and phi=90 degrees.
The solution (there is one!) in this particular case after 4 iterations is F4=2.695233 radians = 154.412 degrees, at which point Y=6.93*10^-11
P.S. I wrote the algorithm in Mathcad to test it, as this is much easier than doing so in Fortran.

[mecej4]

Tony, the first term on the right side of the equation should be sin(phi)*(theta+F-pi) instead of sin(phi)*(theta-F-pi).

[Anthony_Richards]

Yes, that's a misprint in my post - I programmed what you wrote. I will correct my post.

[John_Campbell]

I am interested in the approaches that have been discussed.
Certainly, the use of a quadratic approximation function, such as Mueller's method or a Taylor series expansion will accelerate the convergence, once you are near the solution. I am less convinced of the benefit of Fibonnaci or golden search approaches, as their improvement is probably less effective than a linear approximation approach. (If you analyse these on a statistical approach I think their benefit is outweighed by their complexity.)

My main problem is: do you start from knowing the approximate solution ?  If you look at the (hopefully) attached chart, which graphs Y as a function of F(radians) for Half_angle = 30 deg and theta = 180 deg then if you start for F -ve then your quadratic conversion approach can fail. This is the case when the function is of the form of Sin and you start away from the zero. Anthony's approach of fy(F) = -Y/sin(half_angle) and the derivatives have been plotted to show how these functions vary. If you start near the answer, then a quadratic approximation works well and after a few iterations with most methods you will be near the answer.

What I like about the binary search method is that it is robust and if you start at 2 points where the function is of opposite sign then you should converge to a solution, providing the function is continuous. Once you get close, you could possibly change to a linear or quadratic prediction method, but the problem is deciding when is close enough, especially with complex trig functions.

My first approach at understanding the problem was to graph the function in Excel to see what shape it was then look at where the search would be best or worst applied. It is this estimation of a suitable search range that needs to be improved.
If you start near the answer then any quadratic search rule is going to perform well.

http://www.dropbox.com/s/6rbxkbivnz02kn5/function.png

[DavidWhite]

I haven't attempted to plot the function, so don't know the complications of it.

The search algorithms can be useful but assume that there is only one root in the region searched.  The usefulness comes in where you have a function, for example, which has a minimum just above the axis.  Any method relying on slope in this region will be sent off to extremely large x values, and may end up endlessly hunting for the root.

Search algorthims then become the most efficient method, and as mentioned (partially) before, if the only tolerance of interest is in the x value, then the number of iterations can be determined before starting.  To ensure a tolerance in the y value, the number of iterations can't be determined, however, one could use the linear or quadratic search methods discussed to optimize once in the region of the solution.

[Anthony_Richards]

I think you have to use your knowledge of the equation to tailor your approach. Clearly, the equation presented, containing periodic functions, is almost certain to have multiple roots. Also, my approach would be to limit the magnitude of changes dF0 in F if the computed adjustment is above a certain magnitude. Typically one would apply either an absolute maximum magnitude change, or apply a fixed reduction factor to the computed increment dF0. Also, when no solution exists for the quadratic equation, just apply a maximum value adjustment, with sign taken from the immediately preceding increment or, if the first cycle, choose a sign randomly.

I have introduced such changes to my original algorithm described above and have no trouble with the combination of starting values described by John Campbell. Obviously, any algorithm will be found wanting at some times, and will blow up, but that is why one must utilise knowledge of the equation's basic structure to tweak a basic algorithm to suit. An algorithm that would always converge on every Y=0 point in the 4-dimensional landscape described by (Y,F,theta, phi) would take some pretty exceptional skill which I am nowhere near having.