EDN Admin
Well-known member
I was working on a Solver , and I though that it can be useful to some one else.
So ... here we are.
------------------------------------
Some problems can sometime be relatively easy to solve with a pen and a paper, but when comes time to solve them with a computer, the solution becomes very complicated.
I will start with an example of one of those problem ... Then Ill show a solution.
I choose this example since it is easy to understand (at least the physic of it, if you are not good in math ...)
The problem example is:
Let suppose that you are standing 100 meters above sea level (On a lighthouse may be), looking at the waves right under you.
This is a bad day, and the waves are 3 meters high. And also you notice that a wave hits the lighthouse every 2.1 seconds
Now, suppose that at the moment where a wave is at its lowest, you drop a ball.
The question is: how many waves and fraction of wave will pass before the ball touch the water?
OK, the ball will go down following the equation of an uniformly accelerated object:
F(t) = 100 + (-9.8 * t^2 / 2)
but here the distance that the ball has to travel before touching the water is not constant, it depend on the position of the wave at the time it touch the water.
The height of the wave in a function of the time can be found using the simple harmonic, not damped equation:
F(t) = 3 * cos((2 * pi/2.1) * t + pi)
therefore our solution is given by solving;
100 + (-9.8 * t^2 / 2) = 3 * cos((2 * pi/2.1) * t + pi)
-------------------
The problem here is that it is not possible to isolate the variable "t" out of the equation, so we cannot do
Dim t As Double = [something]
One can always run that in a loop and try to converge to a solution, but if you whish to get a high degree of precision, 6 or 7 digit or more, you will have to do millions of iterations and it will take quite some time before you get to the solution.
-------------------
SO,...
I implemented a class that can solve those system of equations. All you have to do is to make a function for each equation, send that to the solver and it will gives you the answer
Then your code to solve the above problem becomes as simple as this:
Dim Frequency As Double = 2.1
Dim Fnct3 As Func(Of Double, Double) = Function(t As Double) 3 * Math.Cos((Math.PI * 2 / Frequency) * t + Math.PI)
Dim Fnct4 As Func(Of Double, Double) = Function(t As Double) 100 - 9.8 * t ^ 2 / 2
Dim Time() As Double = Solver.Solve(Fnct3, Fnct4)
Dim Nb_Wave As Double = Time(1) / Frequency
And you will find that 2.1671917141288 waves has passed before your ball touch the water
--------------------------
The solver will solve any equation, except:
1) equation of the first degree order F(x) = Ax + B (You probably dont need a solver for that anyway)
2) The function has to be defined at all point of the range. (Obvoiusly, dont try to compute the tangent of pi/2 or something like that ... It never works). We will see about setting the range later
------------------------
The solver is using a quadratic convergence method. Both overloads of the method "Solve" typically return the solutions in less than 0.5 milliseconds. And both overloads of the function SolveFrom returns within 0.1 milliseconds
The solver expose those methods
1) Solve(F As Func(Of Double, Double), G As Func(Of Double, Double)) as double()
This is the function I used in the example above
- It takes as argument the 2 functions in the equation system and return an array of Double containing all the solutions found in the range (We will see later about that range)
2) Solve(F As Func(Of Double, Double)) as double()
- it takes as argument a single function, and return an array of Double containing all the solution found in the range.
- This overload is usefull to solve this type of problem
Find the value of x in F(x) = x*e^(sin(x))
where only one equation has to be solved.
3) SolveFrom(F As Func(Of Double, Double), G As Func(Of Double, Double), From As Double) as double
4) SolveFrom(F As Func(Of Double, Double)), From As Double) as double
- Those two functions, are similar to the previous except that instead of searching for all the solutions, it search for the single solution that is the closest to the value passed in the "From" parameter.
- Those function are usefull in the cases where there are a large number of solutions and we need to select one of them base on a particular value
--------------------
The solver also expose those properties
1) Range: This property contains a structure of 2 values "Left" and "Right"
- Use
A) When searching for a solution, the solver will only return the solution found between those 2 values. (Just think at equation that have an infinite number of solution { ex: Sin(x) = 0 }
B) The range can also be used to restrict the search of solution where the function is defined. { ex: set the range to be between -pi/2+d and pi/2-d if the function is tan(x) }
C) The solver may fail to converge in the unlikely case where "Left+0.0001" or "Right" value of the range would happen to be set where the tangente of the function would be exactly 0.000000000
When this happen, the functions Solve return and empty array and the function SolveFrom return NaN. In this case, it will be required to change the range values and recall the function.
I could have handle this internally in the function, but I choose not, to make the function as quick as possible. And since those problematic values of range are usualy obvious just looking at the function, it is possible to set the range before calling the function.
NB: the functions will also return an empty array / NaN if no solution exists in the given range
NB: By default, the range is set to [-100, 100]
2) ScanningDensity: If you are expecting two solutions or more within a very short interval, or if your range is very large, you may want to increase the value of this property. By default it is set to 100. This should be appropriate for most functions. When this value is set to a value too small, some solutions may be skipped. ( Increasing this value does increase the time the function takes to return)
3) h: this property can be used to fine tune the error of the solutions find by the functions. Typically the error will be less than 0.000000000001, and this property value should not be changed. (In any case between, h=1 to h=0.001 is the logical range of values.)
---------------------
Here is the code for this solver
Friend NotInheritable Class Solver
Public Const Author As String = "Luc Chenier"
Public Shared Property h As Double = 0.1
Public Shared Property Range As FocussedRange = New FocussedRange With {.Left = -100, .Right = 100}
Public Shared Property ScanningDensity As Integer = 100
Public Shared Function Solve(ByVal F As Func(Of Double, Double)) As Double()
Dim Solutions As New List(Of Double)
Dim Left = SolveFrom(F, Range.Left)
If Double.IsNaN(Left) Then Return New Double() {}
Dim right = SolveFrom(F, Range.Right)
If Math.Abs(Left - right) < 0.000000001 Then
Return New Double() {Left}
Else
Dim Neg As Boolean = If(F(Left + 0.0001) < 0, True, False)
Solutions.Add(Left)
For x As Double = Left + 0.0001 To right - 0.0001 Step (right - Left) / ScanningDensity
If (If(F(x) < 0, True, False) <> Neg) Then
Solutions.Add(SolveFrom(F, x))
Neg = Not Neg
End If
Next
Solutions.Add(right)
Return Solutions.ToArray
End If
End Function
Public Shared Function Solve(ByVal F As Func(Of Double, Double), ByVal G As Func(Of Double, Double)) As Double()
Dim Solutions As New List(Of Double)
Dim Left = SolveFrom(F, G, Range.Left)
If Double.IsNaN(Left) Then Return New Double() {}
Dim right = SolveFrom(F, G, Range.Right)
If Math.Abs(Left - right) < 0.000000001 Then
Return New Double() {Left}
Else
Dim Neg As Boolean = If(F(Left + 0.0001) - G(Left + 0.0001) < 0, True, False)
Solutions.Add(Left)
For x As Double = Left + 0.0001 To right - 0.0001 Step (right - Left) / ScanningDensity
If (If((F(Left + 0.0001) - G(Left + 0.0001)) < 0, True, False) <> Neg) Then
Solutions.Add(SolveFrom(F, G, x))
Neg = Not Neg
End If
Next
Solutions.Add(right)
Return Solutions.ToArray
End If
End Function
Public Shared Function SolveFrom(ByVal F As Func(Of Double, Double), ByVal From As Double) As Double
Dim S As Double = From
Dim S1 As Double = 0.0
Dim Iteration As Integer = 0
Do
Iteration += 1
S1 = S - (F(S) / Derivative(F, S))
If Math.Abs(S - S1) < 0.00000000001 Then Exit Do
If Iteration = 100 Then
S = Double.NaN
Exit Do
End If
S = S1
Loop
Return S
End Function
Public Shared Function SolveFrom(ByVal F As Func(Of Double, Double), ByVal G As Func(Of Double, Double), ByVal From As Double) As Double
Dim S As Double = From
Dim S1 As Double = 0.0
Dim Iteration As Integer = 0
Do
Iteration += 1
S1 = S - (F(S) - G(S)) / (Derivative(F, S) - Derivative(G, S))
If Math.Abs(S - S1) < 0.00000000001 Then Exit Do
If Iteration = 100 Then
S = Double.NaN
Exit Do
End If
S = S1
Loop
Return S
End Function
Private Shared Function Derivative(ByVal F As Func(Of Double, Double), ByVal X As Double) As Double
Dim V1 As Double = 3 * F(X - 4 * h)
Dim V2 As Double = -32 * F(X - 3 * h)
Dim V3 As Double = 168 * F(X - 2 * h)
Dim V4 As Double = -672 * F(X - h)
Dim V5 As Double = 672 * F(X + h)
Dim V6 As Double = -168 * F(X + 2 * h)
Dim V7 As Double = 32 * F(X + 3 * h)
Dim V8 As Double = -3 * F(X + 4 * h)
Return (V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8) / (840 * h)
End Function
Structure FocussedRange
Public Left As Double
Public Right As Double
End Structure
End Class
View the full article
So ... here we are.
------------------------------------
Some problems can sometime be relatively easy to solve with a pen and a paper, but when comes time to solve them with a computer, the solution becomes very complicated.
I will start with an example of one of those problem ... Then Ill show a solution.
I choose this example since it is easy to understand (at least the physic of it, if you are not good in math ...)
The problem example is:
Let suppose that you are standing 100 meters above sea level (On a lighthouse may be), looking at the waves right under you.
This is a bad day, and the waves are 3 meters high. And also you notice that a wave hits the lighthouse every 2.1 seconds
Now, suppose that at the moment where a wave is at its lowest, you drop a ball.
The question is: how many waves and fraction of wave will pass before the ball touch the water?
OK, the ball will go down following the equation of an uniformly accelerated object:
F(t) = 100 + (-9.8 * t^2 / 2)
but here the distance that the ball has to travel before touching the water is not constant, it depend on the position of the wave at the time it touch the water.
The height of the wave in a function of the time can be found using the simple harmonic, not damped equation:
F(t) = 3 * cos((2 * pi/2.1) * t + pi)
therefore our solution is given by solving;
100 + (-9.8 * t^2 / 2) = 3 * cos((2 * pi/2.1) * t + pi)
-------------------
The problem here is that it is not possible to isolate the variable "t" out of the equation, so we cannot do
Dim t As Double = [something]
One can always run that in a loop and try to converge to a solution, but if you whish to get a high degree of precision, 6 or 7 digit or more, you will have to do millions of iterations and it will take quite some time before you get to the solution.
-------------------
SO,...
I implemented a class that can solve those system of equations. All you have to do is to make a function for each equation, send that to the solver and it will gives you the answer
Then your code to solve the above problem becomes as simple as this:
Dim Frequency As Double = 2.1
Dim Fnct3 As Func(Of Double, Double) = Function(t As Double) 3 * Math.Cos((Math.PI * 2 / Frequency) * t + Math.PI)
Dim Fnct4 As Func(Of Double, Double) = Function(t As Double) 100 - 9.8 * t ^ 2 / 2
Dim Time() As Double = Solver.Solve(Fnct3, Fnct4)
Dim Nb_Wave As Double = Time(1) / Frequency
And you will find that 2.1671917141288 waves has passed before your ball touch the water
--------------------------
The solver will solve any equation, except:
1) equation of the first degree order F(x) = Ax + B (You probably dont need a solver for that anyway)
2) The function has to be defined at all point of the range. (Obvoiusly, dont try to compute the tangent of pi/2 or something like that ... It never works). We will see about setting the range later
------------------------
The solver is using a quadratic convergence method. Both overloads of the method "Solve" typically return the solutions in less than 0.5 milliseconds. And both overloads of the function SolveFrom returns within 0.1 milliseconds
The solver expose those methods
1) Solve(F As Func(Of Double, Double), G As Func(Of Double, Double)) as double()
This is the function I used in the example above
- It takes as argument the 2 functions in the equation system and return an array of Double containing all the solutions found in the range (We will see later about that range)
2) Solve(F As Func(Of Double, Double)) as double()
- it takes as argument a single function, and return an array of Double containing all the solution found in the range.
- This overload is usefull to solve this type of problem
Find the value of x in F(x) = x*e^(sin(x))
where only one equation has to be solved.
3) SolveFrom(F As Func(Of Double, Double), G As Func(Of Double, Double), From As Double) as double
4) SolveFrom(F As Func(Of Double, Double)), From As Double) as double
- Those two functions, are similar to the previous except that instead of searching for all the solutions, it search for the single solution that is the closest to the value passed in the "From" parameter.
- Those function are usefull in the cases where there are a large number of solutions and we need to select one of them base on a particular value
--------------------
The solver also expose those properties
1) Range: This property contains a structure of 2 values "Left" and "Right"
- Use
A) When searching for a solution, the solver will only return the solution found between those 2 values. (Just think at equation that have an infinite number of solution { ex: Sin(x) = 0 }
B) The range can also be used to restrict the search of solution where the function is defined. { ex: set the range to be between -pi/2+d and pi/2-d if the function is tan(x) }
C) The solver may fail to converge in the unlikely case where "Left+0.0001" or "Right" value of the range would happen to be set where the tangente of the function would be exactly 0.000000000
When this happen, the functions Solve return and empty array and the function SolveFrom return NaN. In this case, it will be required to change the range values and recall the function.
I could have handle this internally in the function, but I choose not, to make the function as quick as possible. And since those problematic values of range are usualy obvious just looking at the function, it is possible to set the range before calling the function.
NB: the functions will also return an empty array / NaN if no solution exists in the given range
NB: By default, the range is set to [-100, 100]
2) ScanningDensity: If you are expecting two solutions or more within a very short interval, or if your range is very large, you may want to increase the value of this property. By default it is set to 100. This should be appropriate for most functions. When this value is set to a value too small, some solutions may be skipped. ( Increasing this value does increase the time the function takes to return)
3) h: this property can be used to fine tune the error of the solutions find by the functions. Typically the error will be less than 0.000000000001, and this property value should not be changed. (In any case between, h=1 to h=0.001 is the logical range of values.)
---------------------
Here is the code for this solver
Friend NotInheritable Class Solver
Public Const Author As String = "Luc Chenier"
Public Shared Property h As Double = 0.1
Public Shared Property Range As FocussedRange = New FocussedRange With {.Left = -100, .Right = 100}
Public Shared Property ScanningDensity As Integer = 100
Public Shared Function Solve(ByVal F As Func(Of Double, Double)) As Double()
Dim Solutions As New List(Of Double)
Dim Left = SolveFrom(F, Range.Left)
If Double.IsNaN(Left) Then Return New Double() {}
Dim right = SolveFrom(F, Range.Right)
If Math.Abs(Left - right) < 0.000000001 Then
Return New Double() {Left}
Else
Dim Neg As Boolean = If(F(Left + 0.0001) < 0, True, False)
Solutions.Add(Left)
For x As Double = Left + 0.0001 To right - 0.0001 Step (right - Left) / ScanningDensity
If (If(F(x) < 0, True, False) <> Neg) Then
Solutions.Add(SolveFrom(F, x))
Neg = Not Neg
End If
Next
Solutions.Add(right)
Return Solutions.ToArray
End If
End Function
Public Shared Function Solve(ByVal F As Func(Of Double, Double), ByVal G As Func(Of Double, Double)) As Double()
Dim Solutions As New List(Of Double)
Dim Left = SolveFrom(F, G, Range.Left)
If Double.IsNaN(Left) Then Return New Double() {}
Dim right = SolveFrom(F, G, Range.Right)
If Math.Abs(Left - right) < 0.000000001 Then
Return New Double() {Left}
Else
Dim Neg As Boolean = If(F(Left + 0.0001) - G(Left + 0.0001) < 0, True, False)
Solutions.Add(Left)
For x As Double = Left + 0.0001 To right - 0.0001 Step (right - Left) / ScanningDensity
If (If((F(Left + 0.0001) - G(Left + 0.0001)) < 0, True, False) <> Neg) Then
Solutions.Add(SolveFrom(F, G, x))
Neg = Not Neg
End If
Next
Solutions.Add(right)
Return Solutions.ToArray
End If
End Function
Public Shared Function SolveFrom(ByVal F As Func(Of Double, Double), ByVal From As Double) As Double
Dim S As Double = From
Dim S1 As Double = 0.0
Dim Iteration As Integer = 0
Do
Iteration += 1
S1 = S - (F(S) / Derivative(F, S))
If Math.Abs(S - S1) < 0.00000000001 Then Exit Do
If Iteration = 100 Then
S = Double.NaN
Exit Do
End If
S = S1
Loop
Return S
End Function
Public Shared Function SolveFrom(ByVal F As Func(Of Double, Double), ByVal G As Func(Of Double, Double), ByVal From As Double) As Double
Dim S As Double = From
Dim S1 As Double = 0.0
Dim Iteration As Integer = 0
Do
Iteration += 1
S1 = S - (F(S) - G(S)) / (Derivative(F, S) - Derivative(G, S))
If Math.Abs(S - S1) < 0.00000000001 Then Exit Do
If Iteration = 100 Then
S = Double.NaN
Exit Do
End If
S = S1
Loop
Return S
End Function
Private Shared Function Derivative(ByVal F As Func(Of Double, Double), ByVal X As Double) As Double
Dim V1 As Double = 3 * F(X - 4 * h)
Dim V2 As Double = -32 * F(X - 3 * h)
Dim V3 As Double = 168 * F(X - 2 * h)
Dim V4 As Double = -672 * F(X - h)
Dim V5 As Double = 672 * F(X + h)
Dim V6 As Double = -168 * F(X + 2 * h)
Dim V7 As Double = 32 * F(X + 3 * h)
Dim V8 As Double = -3 * F(X + 4 * h)
Return (V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8) / (840 * h)
End Function
Structure FocussedRange
Public Left As Double
Public Right As Double
End Structure
End Class
View the full article