Subsections

• Problem:
• Solution:

Ordinary Differential Equations

Problem:

Solve an ordinary differential equation:

y' = f(x, y)

using a numerical (e.g., trapezoidal) method. Is should be understood that there is a mathematical relationship between x and y, e.g., x= g(y), and this relationship is what we want to find out numerically.

Solution:

Numerical methods return functions as an explicit set of pairs (x, y) instead of as an analytic expression. We need to supply the boundaries x₀ and x_n-1 of the interval in which the function is to be calculated, and the number of points n to be taken into account in that interval. Thus, n is the number of integration steps, x₀ and x_n-1 the integration limits, and $h = \frac{x_{n-1} - x_0}{n}$ the integration step--the distance between two consecutive sampling points. Then the equation

$$y_{i+1} = y_i + h \frac{f(x_i, y_i) + f(x_{i+1}, y_{i+1})}{2}$$

can be used to numerically approximate the function y = g(x) (this is the expression of the trapezoidal method of integration).

Note that the expression for y_i+1 involves y_i+1 itself. The effect of the above formula, applied to all the points in the integration segment, is to originate a set of equations; for them to be solved we need at least one y_i to be defined. The top level call gives us access both to y₀ and y_n-1, but setting any y_iwill actually suffice.

Modelling the recurrence equation:

yn1(H, Xi, Yi, Xi1, Yi1):-
   Yi1 = Yi + H*(Fi+Fi1)/2,
   Xi1 = Xi + H,
   f(Xi, Yi, Fi),
   f(Xi1, Yi1, Fi1).

The recurrent loop relates points of y_i = g(x_i) using the previous predicate:

loop(_H, Xn, Yn, Xn, Yn, [Xn], [Yn]).
loop( H, Xi, Yi, Xn, Yn, [Xi|Xs], [Yi|Ys]):-
   lelin(Xi, Xn),
   yn1(H, Xi, Yi, Xi1, Yi1),
   loop(H, Xi1, Yi1, Xn, Yn, Xs, Ys).

In this predicate the first clause implements the stop condition of the loop, and the recursive clause extends the recurrence equation to the selected interval:

• Stop condition: when the current pair (x_i, y_i) equals the boundary condition (x_n, y_n). If the constraint solver is not accurate enough in the mathematical operations, this stop condition could fail: adding $\frac{x_n - x_0}{n}$ with itself n times may not come out with the result x_n - x₀, due to approximation problems. This can be worked around by stopping the recurrence when the current x coordinate has gone beyond the upper boundary.
• The recursive clause relates (x_i, y_i) with (x_i+1, y_i+1) using the recurrence equation.

The top level call computes the integration step, the list of xcoordinates and the corresponding y values, and writes the result. We could return the results as two lists of x coordinates and yvalues, but we will do a little of formatting:

sv(X0, Y0, Xn, Yn, N):-
   H = (Xn - X0)/N,
   loop(H, X0, Y0, Xn, Yn, Xs, Ys),
   write_res(Xs, Ys).

write_res([], []).
write_res([X|Xs], [Y|Ys]):-
   write(y(X) = Y), nl,
   write_res(Xs, Ys).

To solve, for example, y' = -2xy, we just need to define the relationship among x, y and f(x,y):

f(X, Y, -2*X*Y).

(the program above will work for other differential equations just changing this clause). An example call, where the integration bounds are -3 and 3, twenty sampling points where used, and $g(-3) = \frac{1}{2000}$ is the following:

?- sv(-3, 1/2000, 3, Yn, 20).
   y(-3)=1/2000
   y(-27/10)=1/200
   y(-12/5)=181/5600
   y(-21/10)=7783/51800
   y(-9/5)=1268629/2382800
   y(-3/2)=1268629/851000
   y(-6/5)=36790241/10892800
   y(-9/10)=625434097/99396800
   y(-3/5)=79430130319/8150537600
   y(-3/10)=4686377688821/370849460800
   y(0)=510815168081489/37084946080000
   y(3/10)=4686377688821/370849460800
   y(3/5)=79430130319/8150537600
   y(9/10)=625434097/99396800
   y(6/5)=36790241/10892800
   y(3/2)=1268629/851000
   y(9/5)=1268629/2382800
   y(21/10)=7783/51800
   y(12/5)=181/5600
   y(27/10)=1/200
   y(3)=1/2000

   Yn = 1/2000.