# Differential equations

## Description

A simulation program for qualitative differential equations.

Source: PROLOG programming for artificial intelligence, 3rd Edition, Harlow, 2001, ISBN 0-201-40375-7.

Program source code: differential_equations.pl

## Listing

```% Figure 20.8  A simulation program for qualitative differential equations.
% The program uses the usual list predicates member/2 and conc/3.

:- op( 900, fy, not).

% not Goal): negation as failure;
%   Note: This is often available as a built-in predicate,
%   often written as prefix operator "\+", e.g. \+ likes(mary,snakes)

not Goal  :-
Goal, !, fail
;
true.

% An interpreter for Qualitative Differential Equations

:- op( 100, xfx, ..).
:- op( 500, xfx, :).

% landmarks( Domain, [ Land1, Land2, ... ]
%   Land1, Land2 etc. are landmarks for domain Domain
%   This is part of qualitative model definition, user-defined

%  correspond( Constraint):
%    Constraint specifies corresponding values for a type of constr.

correspond( sum( Dom1:zero, Dom2:zero, Dom3:zero)).

correspond( sum( Dom1:L, Dom2:zero, Dom1:L))  :-
qmag( Dom1:L), L \== zero, not (L = _.._).   % L is nonzero landmark in Dom1

correspond( sum( Dom1:zero, Dom2:L, Dom2:L))  :-
qmag( Dom2:L), L \== zero, not (L = _.._).   % L is nonzero landmark in Dom2

correspond( sum( V1, V2, V3))  :-
correspond( V1, V2, V3).                     % User-defined corr. values

% The following is a dummy definition of correspond/3
% just to avoid undefined-predicate error with some Prologs

correspond( dummy, dummy, dummy).

%  qmag( Domain:QualMagnitude)

qmag( Domain:Qm)  :-
landmarks( Domain, Lands),
qmag( Lands, Qm).

qmag( Lands, L)  :-
member( L, Lands),
L \== minf, L \== inf.           % A finite landmark

qmag( Lands, L1..L2)  :-           % Interval
conc( _, [L1,L2 | _], Lands).    % Two adjacent landmarks

%  relative_qmag( Domain1:QM, Domain2:Landmark, Sign):
%    Sign is the sign of the difference between QM and Landmark
%    if QM < Landmark then Sign = neg, etc.

relative_qmag( Domain:Ma..Mb, Domain:Land, Sign)  :-  !,
landmarks( Domain, Lands),
( compare_lands( Ma, Land, Lands, neg), Sign = neg, !
;
Sign = pos
).

relative_qmag( Domain:M1, Domain:M2, Sign)  :-
landmarks( Domain, Lands),
compare_lands( M1, M2, Lands, Sign), !.

% qdir( Qdir, Sign):
%   Qdir is qualitative direction of change with sign Sign

qdir( dec, neg).
qdir( std, zero).
qdir( inc, pos).

% Laws of qualitative summation

%  qsum( Q1, Q2, Q3):
%    Q3 = Q1 + Q2, qualitative sum over domain [pos,zero,neg]

qsum( pos, pos, pos).
qsum( pos, zero, pos).
qsum( pos, neg, pos).
qsum( pos, neg, zero).
qsum( pos, neg, neg).
qsum( zero, pos, pos).
qsum( zero, zero, zero).
qsum( zero, neg, neg).
qsum( neg, pos, pos).
qsum( neg, pos, zero).
qsum( neg, pos, neg).
qsum( neg, zero, neg).
qsum( neg, neg, neg).

%  qdirsum( D1, D2, D3):
%    qualitative sum over directions of change

qdirsum( D1, D2, D3)  :-
qdir( D1, Q1), qdir( D2, Q2), qdir( D3, Q3),
qsum( Q1, Q2, Q3).

% sum( QV1, QV2, QV3):
%   QV1 = QV2 + QV3,
%   qualitative sum over qualitative values of form Domain:Qmag/Dir
%   When called, this predicate assumes that the
%   domains of all three arguments are instantiated

sum( D1:QM1/Dir1, D2:QM2/Dir2, D3:QM3/Dir3)  :-
qdirsum( Dir1, Dir2, Dir3),     % Directions of change: Dir1 + Dir2 = Dir3
qmag( D1:QM1), qmag( D2:QM2), qmag( D3:QM3),
% QM1+QM2=QM3 must be consistent with all corresponding values:
not (
correspond( sum( D1:V1, D2:V2, D3:V3)),     % V1 + V2 = V3
relative_qmag( D1:QM1, D1:V1, Sign1),
relative_qmag( D2:QM2, D2:V2, Sign2),
relative_qmag( D3:QM3, D3:V3, Sign3),
not qsum( Sign1, Sign2, Sign3) ).

%  mplus( X, Y):
%    Y is a monotonically increasing function of X

mplus( D1:QM1/Dir, D2:QM2/Dir)  :-     % Equal directions of change
qmag( D1:QM1), qmag( D2:QM2),
% QM1, QM2 consistent with all corresponding values between D1, D2:
not ( correspond( D1:V1, D2:V2),
relative_qmag( D1:QM1, D1:V1, Sign1),
relative_qmag( D2:QM2, D2:V2, Sign2),
Sign1 \== Sign2 ).

% deriv( Var1, Var2):
%   time derivative of Var1 is qualitatively equal Var2

deriv( Dom1:Qmag1/Dir1, Dom2:Qmag2/Dir2)  :-
qdir( Dir1, Sign1),
qmag( Dom2:Qmag2),
relative_qmag( Dom2:Qmag2, Dom2:zero, Sign2), % Sign2 = sign of Qmag2
Sign1 = Sign2.

% transition( Domain:Qmag1/Dir1, Domain:Qmag2/Dir2):
%   Variable state transitions between "close" time points

transition( Dom:L1..L2/std, Dom:L1..L2/Dir2)  :-
qdir( Dir2, AnySign).

transition( Dom:L1..L2/inc, Dom:L1..L2/inc).

transition( Dom:L1..L2/inc, Dom:L1..L2/std).

transition( Dom:L1..L2/inc, Dom:L2/inc)  :-
L2 \== inf.

transition( Dom:L1..L2/inc, Dom:L2/std)  :-
L2 \== inf.

transition( Dom:L1..L2/dec, Dom:L1..L2/dec).

transition( Dom:L1..L2/dec, Dom:L1..L2/std).

transition( Dom:L1..L2/dec, Dom:L1/dec)  :-
L1 \== minf.

transition( Dom:L1..L2/dec, Dom:L1/std)  :-
L1 \== minf.

transition( Dom:L1/std, Dom:L1/std)  :-
L1 \== A..B.    % L1 not an interval

transition( Dom:L1/std, Dom:L1..L2/inc)  :-
qmag( Dom:L1..L2).

transition( Dom:L1/std, Dom:L0..L1/dec)  :-
qmag( Dom:L0..L1).

transition( Dom:L1/inc, Dom:L1..L2/inc)  :-
qmag( Dom:L1..L2).

transition( Dom:L1/dec, Dom:L0..L1/dec)  :-
qmag( Dom:L0..L1).

% system_trans( State1, State2):
%   System state transition;
%   system state is a list of variable values

system_trans( [], []).

system_trans( [Val1 | Vals1], [Val2 | Vals2])  :-
transition( Val1, Val2),
system_trans( Vals1, Vals2).

% legal_trans( State1, State2):
%   possible transition between states according to model

legal_trans( State1, State2)  :-
system_trans( State1, State2),
State1 \== State2,            % Qualitatively different next state
legalstate( State2).          % Legal according to model

% simulate( SystemStates, MaxLength):
%  SystemStates is a sequence of states of simulated system
%  not longer than MaxLength

simulate( [State], MaxLength)  :-
( MaxLength = 1                % Max length reached
;
not legal_trans( State, _)   % No legal next state
) , !.

simulate( [State1,State2 | Rest], MaxLength)  :-
MaxLength > 1, NewMaxL is MaxLength - 1,
legal_trans( State1, State2),
simulate( [State2 | Rest], NewMaxL).

% simulate( InitialState, QualBehaviour, MaxLength)

simulate( InitialState, [InitialState | Rest], MaxLength)  :-
legalstate( InitialState),             % Satisfy system's model
simulate( [InitialState | Rest], MaxLength).

% compare_lands( X1, X2, List, Sign):
%   if X1 before X2 in List then Sign = neg
%   if X2 before X1 then Sign = pos else Sign = zero

compare_lands( X1, X2, [First | Rest], Sign)  :-
X1 = X2, !, Sign = zero
;
X1 = First, !, Sign = neg
;
X2 = First, !, Sign = pos
;
compare_lands( X1, X2, Rest, Sign).

% member( X, List): X is a member of List

member(X,[X | _]).                 % X is head of list

member( X, [_ | Rest])  :-
member( X, Rest).                % X is in body of list

%  conc(L1,L2,L3): list L3 is th econcatenation of lists L1 and L2

conc( [], L, L).

conc( [X | L1], L2, [X | L3])  :-
conc( L1, L2, L3).```