42aabce1bb
git-svn-id: https://yap.svn.sf.net/svnroot/yap/trunk@1973 b08c6af1-5177-4d33-ba66-4b1c6b8b522a
264 lines
5.7 KiB
Plaintext
264 lines
5.7 KiB
Plaintext
|
|
:- protocol(find_rootp).
|
|
|
|
:- info([
|
|
version is 1.1,
|
|
author is 'Paulo Moura and Paulo Nunes',
|
|
date is 2006/11/26,
|
|
comment is 'Default protocol for root find algorithms.']).
|
|
|
|
:- public(find_root/5).
|
|
:- mode(find_root(+object_identifier, +float, +float, +float, -float), one).
|
|
:- info(find_root/5, [
|
|
comment is 'Find the root of a function in the interval [A, B] given a maximum aproximation error.',
|
|
argnames is ['Function', 'A', 'B', 'Error', 'Zero']]).
|
|
|
|
:- end_protocol.
|
|
|
|
|
|
|
|
:- protocol(functionp).
|
|
|
|
:- info([
|
|
version is 1.1,
|
|
author is 'Paulo Moura and Paulo Nunes',
|
|
date is 2006/11/26,
|
|
comment is 'Default protocol for real functions of a single real variable.']).
|
|
|
|
:- public(eval/2).
|
|
:- mode(eval(+float, -float), one).
|
|
:- info(eval/2, [
|
|
comment is 'Calculates the function value.',
|
|
argnames is ['X', 'Fx']]).
|
|
|
|
:- public(evald/2).
|
|
:- mode(evald(+float, -float), one).
|
|
:- info(evald/2, [
|
|
comment is 'Calculates the value of the function derivative.',
|
|
argnames is ['X', 'DFx']]).
|
|
|
|
:- end_protocol.
|
|
|
|
|
|
|
|
:- object(f1,
|
|
implements(functionp)).
|
|
|
|
% x^2 - 4
|
|
% 2.0
|
|
|
|
eval(X, Y) :-
|
|
Y is X * X - 4.
|
|
|
|
evald(X, Y) :-
|
|
Y is 2 * X.
|
|
|
|
:- end_object.
|
|
|
|
|
|
|
|
:- object(f2,
|
|
implements(functionp)).
|
|
|
|
% x^7 + 9x^5 - 13x - 17
|
|
% 1.29999999999945448
|
|
|
|
eval(X, Y) :-
|
|
Y is X**7 + 9*X**5 - 13*X - 17.
|
|
|
|
evald(X, Y) :-
|
|
Y is 7*X**6 + 45*X**4 - 13.
|
|
|
|
:- end_object.
|
|
|
|
|
|
|
|
:- object(f3,
|
|
implements(functionp)).
|
|
|
|
% (x - sqrt(2))^7
|
|
% 1.41421356237309537
|
|
|
|
eval(X, Y) :-
|
|
Y is (X - sqrt(2.0))**8.
|
|
|
|
evald(X, Y) :-
|
|
Y is 8*(X - sqrt(2.0))**7.
|
|
|
|
:- end_object.
|
|
|
|
|
|
|
|
:- object(f4,
|
|
implements(functionp)).
|
|
|
|
% x + x^2*sin(2.0/x)
|
|
% 0.0
|
|
|
|
eval(X, Y) :-
|
|
Y is X + (X**2)*sin(2.0/X).
|
|
|
|
evald(X, Y) :-
|
|
Y is 1 + 2*X*sin(2.0/X) - 2*cos(2.0/X).
|
|
|
|
:- end_object.
|
|
|
|
|
|
|
|
:- object(bisection,
|
|
implements(find_rootp)).
|
|
|
|
:- info([
|
|
version is 1.2,
|
|
author is 'Paulo Moura and Paulo Nunes',
|
|
date is 2007/7/7,
|
|
comment is 'Bisection algorithm.']).
|
|
|
|
find_root(Function, A, B, Error, Zero) :-
|
|
Function::eval(A, Fa),
|
|
Function::eval(B, Fb),
|
|
( Fa > 0.0, Fb < 0.0 ->
|
|
true
|
|
; Fa < 0.0, Fb > 0.0
|
|
),
|
|
X0 is (A + B) / 2.0,
|
|
Function::eval(X0, F0),
|
|
bisection(Function, A, B, X0, F0, Error, Zero).
|
|
|
|
bisection(_, _, _, Xn, Fn, Error, Xn) :-
|
|
abs(Fn) < Error,
|
|
!.
|
|
bisection(Function, An, Bn, _, _, Error, Zero) :-
|
|
Xn1 is (An + Bn) / 2.0,
|
|
Function::eval(Xn1, Fn1),
|
|
Function::eval(An, FAn),
|
|
( Fn1*FAn < 0.0 ->
|
|
An1 is An,
|
|
Bn1 is Xn1
|
|
; An1 is Xn1,
|
|
Bn1 is Bn
|
|
),
|
|
bisection(Function, An1, Bn1, Xn1, Fn1, Error, Zero).
|
|
|
|
:- end_object.
|
|
|
|
|
|
|
|
:- object(newton,
|
|
implements(find_rootp)).
|
|
|
|
:- info([
|
|
version is 1.2,
|
|
author is 'Paul Crocker... No More Coffee',
|
|
date is 2007/07/06,
|
|
comment is 'Newton algorithm.']).
|
|
|
|
find_root(Function, Xa, Xb, Deviation, Zero) :-
|
|
Ac is (Xb - Xa) / 2,
|
|
newton(Function, Xa, Ac, Deviation, Zero).
|
|
|
|
newton(_, Zero, Ac, Deviation, Zero) :-
|
|
abs(Ac) < Deviation,
|
|
!.
|
|
newton(Function, X0, Ac, Deviation, Zero):-
|
|
Xn1 is X0 + Ac,
|
|
Function::eval(Xn1, Fx),
|
|
Function::evald(Xn1, DFx),
|
|
Ac2 is -(Fx/DFx),
|
|
newton(Function, Xn1, Ac2, Deviation, Zero).
|
|
|
|
:- end_object.
|
|
|
|
|
|
:- object(muller,
|
|
implements(find_rootp)).
|
|
|
|
:- info([
|
|
version is 1.2,
|
|
author is 'Paulo Moura and Paulo Nunes',
|
|
date is 2006/11/26,
|
|
comment is 'Muller algorithm.']).
|
|
|
|
find_root(Function, Xa, Xb, Deviation, Zero) :-
|
|
Xc is (Xa + Xb) / 2.0,
|
|
muller(Function, Xa, Xc, Xb, Deviation, Zero).
|
|
|
|
muller(Function, Xa, Xb, Xc, Deviation, Zero) :-
|
|
Function::eval(Xa, Ya),
|
|
Function::eval(Xb, Yb),
|
|
Function::eval(Xc, Yc),
|
|
H1 is Xb - Xa,
|
|
DDba is (Yb - Ya) / H1,
|
|
Ac is Deviation + 1.0,
|
|
muller(Function, Xa, Xb, Xc, Deviation, Ya, Yb, Yc, Ac, H1, DDba, Zero).
|
|
|
|
muller(_, _, _, Xc, Deviation, _, _, _, Ac, _, _, Xc) :-
|
|
abs(Ac) < Deviation,
|
|
!.
|
|
muller(Function, Xa, Xb, Xc, Deviation, _, Yb, Yc, _, _, DDba, Zero) :-
|
|
H2n is Xc - Xb,
|
|
DDcbn is (Yc - Yb) / H2n,
|
|
Cn is (DDcbn - DDba) / (Xc - Xa),
|
|
Bn is DDcbn + H2n * Cn,
|
|
Rn is Bn * Bn - 4.0 * Yc * Cn,
|
|
( Rn < 0.0 ->
|
|
fail
|
|
; V is sqrt(Rn)
|
|
),
|
|
( Bn > 0.0 ->
|
|
Dn is Bn + V
|
|
; Dn is Bn - V
|
|
),
|
|
Acn is -(2 * Yc / Dn),
|
|
Xan is Xb,
|
|
Xbn is Xc,
|
|
Xcn is Xc + Acn,
|
|
Yan is Yb,
|
|
Ybn is Yc,
|
|
Function::eval(Xcn, Ycn),
|
|
H1n is H2n,
|
|
DDban is DDcbn,
|
|
muller(Function, Xan, Xbn, Xcn, Deviation, Yan, Ybn, Ycn, Acn, H1n, DDban, Zero).
|
|
|
|
:- end_object.
|
|
|
|
|
|
|
|
:- object(function_root,
|
|
implements(find_rootp)).
|
|
|
|
:- info([
|
|
version is 2.0,
|
|
author is 'Paulo Moura and Paulo Nunes',
|
|
date is 2007/07/05,
|
|
comment is 'Multi-threading interface to root finding algorithms.']).
|
|
|
|
:- threaded.
|
|
|
|
:- public(find_root/6).
|
|
:- mode(find_root(+object_identifier, +float, +float, +float, -float, -object_identifier), one).
|
|
:- info(find_root/6, [
|
|
comment is 'Finds the root of a function in the interval [A, B] given a maximum aproximation error. Returns the method used.',
|
|
argnames is ['Function', 'A', 'B', 'Error', 'Zero', 'Method']]).
|
|
|
|
find_root(Function, A, B, Error, Zero, Algorithm) :-
|
|
threaded((
|
|
(bisection::find_root(Function, A, B, Error, Zero), Algorithm = bisection)
|
|
; (newton::find_root(Function, A, B, Error, Zero), Algorithm = newton)
|
|
; (muller::find_root(Function, A, B, Error, Zero), Algorithm = muller)
|
|
)).
|
|
|
|
% an alternative, possibly better definition would be to make the methods simply fail in case of error:
|
|
%
|
|
% find_root(Function, A, B, Error, Zero, Algorithm) :-
|
|
% threaded((
|
|
% (catch(bisection::find_root(Function, A, B, Error, Zero), _, fail), Algorithm = bisection)
|
|
% ; (catch(newton::find_root(Function, A, B, Error, Zero), _, fail), Algorithm = newton)
|
|
% ; (catch(muller::find_root(Function, A, B, Error, Zero), _, fail), Algorithm = muller)
|
|
% )).
|
|
|
|
find_root(Function, A, B, Error, Zero) :-
|
|
find_root(Function, A, B, Error, Zero, _).
|
|
|
|
:- end_object.
|