line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
# ########################################################################################
|
2
|
|
|
|
|
|
|
# A NEWTON RAPHSON OBJECT
|
3
|
|
|
|
|
|
|
# This module takes an equation in symbolic form and uses the Newton Raphson technique
|
4
|
|
|
|
|
|
|
# to solve it.
|
5
|
|
|
|
|
|
|
# Copyright (C) Jonathan Worthington 2004
|
6
|
|
|
|
|
|
|
# This module may be used and distributed under the same terms as Perl.
|
7
|
|
|
|
|
|
|
# ########################################################################################
|
8
|
|
|
|
|
|
|
|
9
|
|
|
|
|
|
|
package Math::Calculus::NewtonRaphson;
|
10
|
1
|
|
|
1
|
|
8565
|
use Math::Calculus::Expression;
|
|
0
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
11
|
|
|
|
|
|
|
use Math::Calculus::Differentiate;
|
12
|
|
|
|
|
|
|
use strict;
|
13
|
|
|
|
|
|
|
our $VERSION = '0.1';
|
14
|
|
|
|
|
|
|
our @ISA = qw/Math::Calculus::Expression/;
|
15
|
|
|
|
|
|
|
our $MAXITERATIONS = 100;
|
16
|
|
|
|
|
|
|
|
17
|
|
|
|
|
|
|
=head1 NAME
|
18
|
|
|
|
|
|
|
|
19
|
|
|
|
|
|
|
Math::Calculus::NewtonRaphson - Algebraic Newton Raphson Implementation
|
20
|
|
|
|
|
|
|
|
21
|
|
|
|
|
|
|
=head1 SYNOPSIS
|
22
|
|
|
|
|
|
|
|
23
|
|
|
|
|
|
|
use Math::Calculus::NewtonRaphson;
|
24
|
|
|
|
|
|
|
|
25
|
|
|
|
|
|
|
# Create an object.
|
26
|
|
|
|
|
|
|
my $exp = Math::Calculus::NewtonRaphson->new;
|
27
|
|
|
|
|
|
|
|
28
|
|
|
|
|
|
|
# Set a variable and expression.
|
29
|
|
|
|
|
|
|
$exp->addVariable('x');
|
30
|
|
|
|
|
|
|
$exp->setExpression('x^2 - 5') or die $exp->getError;
|
31
|
|
|
|
|
|
|
|
32
|
|
|
|
|
|
|
# Apply Newton Raphson.
|
33
|
|
|
|
|
|
|
my $result = $exp->newtonRaphson(2) or die $exp->getError;
|
34
|
|
|
|
|
|
|
print $result; # Prints 1.4142...
|
35
|
|
|
|
|
|
|
|
36
|
|
|
|
|
|
|
|
37
|
|
|
|
|
|
|
=head1 DESCRIPTION
|
38
|
|
|
|
|
|
|
|
39
|
|
|
|
|
|
|
This module can take an algebraic expression, parses it and then uses the Newton Raphson
|
40
|
|
|
|
|
|
|
method to solve the it. The Newton Raphson method relies on the fact that the expression
|
41
|
|
|
|
|
|
|
you pass in evaluates to zero where there is a solution. That is, to solve:-
|
42
|
|
|
|
|
|
|
|
43
|
|
|
|
|
|
|
x^2 = 5
|
44
|
|
|
|
|
|
|
|
45
|
|
|
|
|
|
|
You would need to pass in:-
|
46
|
|
|
|
|
|
|
|
47
|
|
|
|
|
|
|
x^2 - 5
|
48
|
|
|
|
|
|
|
|
49
|
|
|
|
|
|
|
It understands expressions containing any of the operators +, -, *, / and ^ (raise to
|
50
|
|
|
|
|
|
|
power), bracketed expressions to enable correct precedence and the functions ln,
|
51
|
|
|
|
|
|
|
exp, sin, cos, tan, sec, cosec, cot, sinh, cosh, tanh, sech, cosech, coth, asin,
|
52
|
|
|
|
|
|
|
acos, atan, asinh, acosh and atanh.
|
53
|
|
|
|
|
|
|
|
54
|
|
|
|
|
|
|
=head1 EXPORT
|
55
|
|
|
|
|
|
|
|
56
|
|
|
|
|
|
|
None by default.
|
57
|
|
|
|
|
|
|
|
58
|
|
|
|
|
|
|
=head1 METHODS
|
59
|
|
|
|
|
|
|
|
60
|
|
|
|
|
|
|
=cut
|
61
|
|
|
|
|
|
|
|
62
|
|
|
|
|
|
|
# Constructor
|
63
|
|
|
|
|
|
|
# ###########
|
64
|
|
|
|
|
|
|
|
65
|
|
|
|
|
|
|
=item new
|
66
|
|
|
|
|
|
|
|
67
|
|
|
|
|
|
|
$exp = Math::Calculus::NewtonRaphson->new;
|
68
|
|
|
|
|
|
|
|
69
|
|
|
|
|
|
|
Creates a new instance of the Newton Raphson object, which can hold an individual
|
70
|
|
|
|
|
|
|
expression.
|
71
|
|
|
|
|
|
|
|
72
|
|
|
|
|
|
|
=item addVariable
|
73
|
|
|
|
|
|
|
|
74
|
|
|
|
|
|
|
$exp->addVariable('x');
|
75
|
|
|
|
|
|
|
|
76
|
|
|
|
|
|
|
Sets a certain named value in the expression as being a variable. A named value must be
|
77
|
|
|
|
|
|
|
an alphabetic chracter.
|
78
|
|
|
|
|
|
|
|
79
|
|
|
|
|
|
|
=item setExpression
|
80
|
|
|
|
|
|
|
|
81
|
|
|
|
|
|
|
$exp->setExpression('x^2 + 5*x);
|
82
|
|
|
|
|
|
|
|
83
|
|
|
|
|
|
|
Takes an expression in human-readable form and stores it internally as a tree structure,
|
84
|
|
|
|
|
|
|
checking it is a valid expression that the module can understand in the process. Note that
|
85
|
|
|
|
|
|
|
the engine is strict about syntax. For example, note above that you must write 5*x and not
|
86
|
|
|
|
|
|
|
just 5x. Whitespace is allowed in the expression, but does not have any effect on precedence.
|
87
|
|
|
|
|
|
|
If you require control of precedence, use brackets; bracketed expressions will always be
|
88
|
|
|
|
|
|
|
evaluated first, as you would normally expect. The module follows the BODMAS precedence
|
89
|
|
|
|
|
|
|
convention. Returns undef on failure and a true value on success.
|
90
|
|
|
|
|
|
|
|
91
|
|
|
|
|
|
|
=item getExpression
|
92
|
|
|
|
|
|
|
|
93
|
|
|
|
|
|
|
$expr = $exp->getExpression;
|
94
|
|
|
|
|
|
|
|
95
|
|
|
|
|
|
|
Returns a textaul, human readable representation of the expression that is being stored.
|
96
|
|
|
|
|
|
|
|
97
|
|
|
|
|
|
|
=cut
|
98
|
|
|
|
|
|
|
|
99
|
|
|
|
|
|
|
|
100
|
|
|
|
|
|
|
# Newton Raphson.
|
101
|
|
|
|
|
|
|
# ###############
|
102
|
|
|
|
|
|
|
|
103
|
|
|
|
|
|
|
=item newtonRaphson
|
104
|
|
|
|
|
|
|
|
105
|
|
|
|
|
|
|
$result = $exp->newtonRaphson($variable, $guess, %mappings);
|
106
|
|
|
|
|
|
|
|
107
|
|
|
|
|
|
|
Attempts to solve the expression for the given variable using the Newton Raphson method,
|
108
|
|
|
|
|
|
|
using the passed value as the first guess. The mappings hash allows any other non-numeric
|
109
|
|
|
|
|
|
|
constants to be mapped to numeric values - a pre-requisite for solving such equations.
|
110
|
|
|
|
|
|
|
|
111
|
|
|
|
|
|
|
=cut
|
112
|
|
|
|
|
|
|
|
113
|
|
|
|
|
|
|
sub newtonRaphson {
|
114
|
|
|
|
|
|
|
# Get invocant and variable.
|
115
|
|
|
|
|
|
|
my ($self, $variable, $guess, %mappings) = @_;
|
116
|
|
|
|
|
|
|
|
117
|
|
|
|
|
|
|
# Clear error and traceback.
|
118
|
|
|
|
|
|
|
$self->{'error'} = $self->{'traceback'} = '';
|
119
|
|
|
|
|
|
|
|
120
|
|
|
|
|
|
|
# Check variable is in the list of variables.
|
121
|
|
|
|
|
|
|
unless (grep { $_ eq $variable } @{$self->{'variables'}})
|
122
|
|
|
|
|
|
|
{
|
123
|
|
|
|
|
|
|
$self->{'error'} = 'Function variable was not declared.';
|
124
|
|
|
|
|
|
|
return undef;
|
125
|
|
|
|
|
|
|
}
|
126
|
|
|
|
|
|
|
|
127
|
|
|
|
|
|
|
# Attempt to differentiate the expression.
|
128
|
|
|
|
|
|
|
my $diffExp = Math::Calculus::Differentiate->new;
|
129
|
|
|
|
|
|
|
$diffExp->setExpression($self->getExpression);
|
130
|
|
|
|
|
|
|
$diffExp->addVariable($_) foreach @{$self->{'variables'}};
|
131
|
|
|
|
|
|
|
unless ($diffExp->differentiate($variable)) {
|
132
|
|
|
|
|
|
|
$self->{'error'} = 'Unable to differentiate expression';
|
133
|
|
|
|
|
|
|
return undef;
|
134
|
|
|
|
|
|
|
}
|
135
|
|
|
|
|
|
|
|
136
|
|
|
|
|
|
|
# Build up an expression for us to plug values into.
|
137
|
|
|
|
|
|
|
my $fiter = {
|
138
|
|
|
|
|
|
|
operation => '/',
|
139
|
|
|
|
|
|
|
operand1 => $self->{'expression'},
|
140
|
|
|
|
|
|
|
operand2 => $diffExp->getExpressionTree
|
141
|
|
|
|
|
|
|
};
|
142
|
|
|
|
|
|
|
|
143
|
|
|
|
|
|
|
# Now iterate.
|
144
|
|
|
|
|
|
|
my $curGuess = $guess;
|
145
|
|
|
|
|
|
|
my $lastGuess = !$guess;
|
146
|
|
|
|
|
|
|
my $iterations = 0;
|
147
|
|
|
|
|
|
|
while ($iterations < $MAXITERATIONS && $curGuess != $lastGuess) {
|
148
|
|
|
|
|
|
|
# Write traceback.
|
149
|
|
|
|
|
|
|
$self->{'traceback'} .= "$iterations\t$curGuess\n";
|
150
|
|
|
|
|
|
|
|
151
|
|
|
|
|
|
|
# Sub value in.
|
152
|
|
|
|
|
|
|
$lastGuess = $curGuess;
|
153
|
|
|
|
|
|
|
eval {
|
154
|
|
|
|
|
|
|
$curGuess = $lastGuess - $self->evaluateTree($fiter, $variable => $lastGuess, %mappings);
|
155
|
|
|
|
|
|
|
} || ($self->{'error'} ||= "Fatal error! $@");
|
156
|
|
|
|
|
|
|
|
157
|
|
|
|
|
|
|
# Increment iterations counter.
|
158
|
|
|
|
|
|
|
$iterations++;
|
159
|
|
|
|
|
|
|
}
|
160
|
|
|
|
|
|
|
|
161
|
|
|
|
|
|
|
# Return an appropriate value (or lack thereof...).
|
162
|
|
|
|
|
|
|
if ($self->{'error'}) {
|
163
|
|
|
|
|
|
|
return undef;
|
164
|
|
|
|
|
|
|
} else {
|
165
|
|
|
|
|
|
|
return $curGuess;
|
166
|
|
|
|
|
|
|
}
|
167
|
|
|
|
|
|
|
}
|
168
|
|
|
|
|
|
|
|
169
|
|
|
|
|
|
|
|
170
|
|
|
|
|
|
|
=item getTraceback
|
171
|
|
|
|
|
|
|
|
172
|
|
|
|
|
|
|
$exp->getTraceback;
|
173
|
|
|
|
|
|
|
|
174
|
|
|
|
|
|
|
When setExpression and newtonRaphson are called, a traceback is generated to describe
|
175
|
|
|
|
|
|
|
what these functions did. If an error occurs, this traceback can be extremely useful
|
176
|
|
|
|
|
|
|
in helping track down the source of the error.
|
177
|
|
|
|
|
|
|
|
178
|
|
|
|
|
|
|
=item getError
|
179
|
|
|
|
|
|
|
|
180
|
|
|
|
|
|
|
$exp->getError;
|
181
|
|
|
|
|
|
|
|
182
|
|
|
|
|
|
|
When any method other than getTraceback is called, the error message stored is cleared, and
|
183
|
|
|
|
|
|
|
then any errors that occur during the execution of the method are stored. If failure occurs,
|
184
|
|
|
|
|
|
|
call this method to get a textual representation of the error.
|
185
|
|
|
|
|
|
|
|
186
|
|
|
|
|
|
|
=head1 SEE ALSO
|
187
|
|
|
|
|
|
|
|
188
|
|
|
|
|
|
|
The author of this module has a website at L, which has
|
189
|
|
|
|
|
|
|
the latest news about the module and a web-based frontend to allow you to test the module
|
190
|
|
|
|
|
|
|
out for yourself.
|
191
|
|
|
|
|
|
|
|
192
|
|
|
|
|
|
|
=head1 AUTHOR
|
193
|
|
|
|
|
|
|
|
194
|
|
|
|
|
|
|
Jonathan Worthington, Ejonathan@jwcs.netE
|
195
|
|
|
|
|
|
|
|
196
|
|
|
|
|
|
|
=head1 COPYRIGHT AND LICENSE
|
197
|
|
|
|
|
|
|
|
198
|
|
|
|
|
|
|
Copyright (C) 2004 by Jonathan Worthington
|
199
|
|
|
|
|
|
|
|
200
|
|
|
|
|
|
|
This library is free software; you can redistribute it and/or modify
|
201
|
|
|
|
|
|
|
it under the same terms as Perl itself, either Perl version 5.8.1 or,
|
202
|
|
|
|
|
|
|
at your option, any later version of Perl 5 you may have available.
|
203
|
|
|
|
|
|
|
|
204
|
|
|
|
|
|
|
=cut
|
205
|
|
|
|
|
|
|
|
206
|
|
|
|
|
|
|
1;
|
207
|
|
|
|
|
|
|
|