001/* 002 * Shredzone Commons - suncalc 003 * 004 * Copyright (C) 2018 Richard "Shred" Körber 005 * http://commons.shredzone.org 006 * 007 * Licensed under the Apache License, Version 2.0 (the "License"); 008 * you may not use this file except in compliance with the License. 009 * 010 * This program is distributed in the hope that it will be useful, 011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 013 */ 014package org.shredzone.commons.suncalc.util; 015 016import static java.lang.Math.abs; 017 018/** 019 * Finds the root of a function by using the Pegasus method. 020 * 021 * @see <a href="https://en.wikipedia.org/wiki/False_position_method">regula falsi</a> 022 * @since 2.3 023 */ 024public class Pegasus { 025 026 private static final int MAX_ITERATIONS = 30; 027 028 /** 029 * Find the root of the given function within the boundaries. 030 * 031 * @param lower 032 * Lower boundary 033 * @param upper 034 * Upper boundary 035 * @param accuracy 036 * Desired accuracy 037 * @param f 038 * Function to be used for calculation 039 * @return root that was found 040 * @throws ArithmeticException 041 * if the root could not be found in the given accuracy within 042 * {@value #MAX_ITERATIONS} iterations. 043 */ 044 public static Double calculate(double lower, double upper, double accuracy, Function f) { 045 double x1 = lower; 046 double x2 = upper; 047 048 double f1 = f.apply(x1); 049 double f2 = f.apply(x2); 050 051 if (f1 * f2 >= 0.0) { 052 throw new ArithmeticException("No root within the given boundaries"); 053 } 054 055 int i = MAX_ITERATIONS; 056 057 while (i-- > 0) { 058 double x3 = x2 - f2 / ((f2 - f1) / (x2 - x1)); 059 double f3 = f.apply(x3); 060 061 if (f3 * f2 <= 0.0) { 062 x1 = x2; 063 f1 = f2; 064 x2 = x3; 065 f2 = f3; 066 } else { 067 f1 = f1 * f2 / (f2 + f3); 068 x2 = x3; 069 f2 = f3; 070 } 071 072 if (abs(x2 - x1) <= accuracy) { 073 return abs(f1) < abs(f2) ? x1 : x2; 074 } 075 } 076 077 throw new ArithmeticException("Maximum number of iterations exceeded"); 078 } 079 080 /** 081 * The function that is to be solved. 082 * 083 * @since 2.3 084 */ 085 public interface Function { 086 087 /** 088 * Calculate the function result for x. 089 * 090 * @param x 091 * x 092 * @return f(x) 093 */ 094 double apply(double x); 095 096 } 097 098}