# Estimating the parameters of a system of differential equations from inaccurate observations

## remark

This work is based on material licensed under the Creative Commons Attribution 4.0 (CC BY 4.0) International License. Source material: “ESTIMATION OF THE PARAMETERS OF DIFFERENTIAL EQUATIONS FROM INEXACT OBSERVATIONS” by G. Sh. Tsitsiashvili and M. A. Osipov, available at https://journals.vsu.ru/sait/article/view/10247/9468 (date of access: 07/14/2023). No changes were made to the original algorithm/method, but only its implementation was presented and an experiment was carried out.

I note that this work (code) uses a string calculator that I previously created, the operation of which was described in a previous article (https://habr.com/ru/articles/747920/). In this regard, we will not dwell on its functions in this article.

## Motivation and about the article

### The article mentioned above describes a method that allows one to estimate the parameters of a system of differential equations. Since this article aroused my interest, I decided to conduct a numerical experiment in the Scala language, which is closest to me, in order to see real results and get to know the algorithm in more detail.

This is what the Lorentz system looks like, consisting of three differential equations. Usually it is described in such notation, but I will rewrite it so that the notation matches those given in the article and displayed in the code.

Here `yi`

represent variables, and `bi`

are the constants we are going to evaluate from the observations.

The article considers such systems, and it is easy to verify that it corresponds to the Lorentz system (m = 3).

It is stated that if there are some inaccurate observations yi’ with a time step hi,

That

From the system, *10, estimates of our parameters bi are derived, which are approximately equal to bi0

## The first step is data preparation.

The first thing I will do is define the variables we are setting and define the data formats.

```
object LorenzAttractor extends App {
val n = 1000 // *9
val alpha = 1.25 // *8
val h = Math.pow(n, -alpha) // *8
val bi: List[Double] = List(2, 3, 1) // фиксированные параметры bi
val yi0: List[Double] = List(1, 2, 1) // начальне значения во временом наге равном 0
val systems: List[String] = List("b1*(y2-y1)", "y1*(b2-y3)-y2", "y1*y2-b3*y3") // сисиема уравнений *2
val reverseSystems: List[String] = List("F1/(y2-y1)", "(F2+y2)/y1+y3", "(y1*y2 - F3)/y3") // решение системы *10
}
```

`reverseSystem`

– this is a necessary crutch, since I have not yet figured out how to solve arbitrary systems of equations (perhaps my future article will be devoted to this issue). Matrix solutions are suitable only for systems of linear equations.

For further work, it is also necessary to set up a string calculator so that it can work with all the operations that are in `systems`

And `reverseSystems`

namely: “+”, “-“, “*”, “/”.

```
object LorenzAttractor extends App with StringCalculator with ConstantFabric with BinaryOperationFabric {
addBinaryOperation(createBinary("+", 1, (left, right) => left + right, 0))
addBinaryOperation(createBinary("-", 1, (left, right) => left - right, 0))
addBinaryOperation(createBinary("*", 2, (left, right) => left * right, 1))
addBinaryOperation(createBinary("/", 2, (left, right) => left / right, 1))
...
bi.zipWithIndex.foreach(b_index => addConstant(createConstant(s"b${b_index._2 + 1}", b_index._1)))
}
```

For this we use classes `StringCalculator`

, `ConstantFabric`

And`BinaryOperationFabric`

(see my article about the string calculator) and add the appropriate operations.

To add `bi`

into equations, we represent them as constants in a string calculator.

## Step 2 – generating inaccurate observations

first, let’s make a small method that will calculate the values of functions with our variables yk, i

```
private def calkFun(fun: String, values: List[Double]): Double = {
calculate(values.indices.foldLeft(fun)((acc, i) => acc.replace(s"y${i + 1}", s"(${values(i)})")))
}
```

This method replaces all `yi`

in formulas for the corresponding values from `values`

and then calculates the value of the function at the point `Mi(yk,1, ..., yk,m)`

(method `calculate`

– part of the string calculator).

Since the Lorentz *1,2 system does not have an analytical solution, we will solve the system numerically. This is beneficial in a sense, since we need Fuzzy Observations.

As a numerical method, we will use the most elementary method – the Euler method. To understand exactly how I will use it, one formula is enough.

Let us separately generate observations for positive `k`

.

```
val positivePreciseObservations: List[List[Double]] = (1 to n).foldLeft(List(yi0))((acc, itr) => {
acc :+ systems.zipWithIndex.map(fun_index => acc.last(fun_index._2) + h * calkFun(fun_index._1, acc.last))
})
```

The algorithm of this code is as follows:

A variable is created

`positivePreciseObservations`

which is initialized with a list containing the initial observation`yi0`

as the first element. This list will contain a sequence of precise observations of the system.Then the method is applied

`foldLeft`

to the range from 1 to`n`

. In every iteration`itr`

represents the current iteration number.Within each iteration, the following happens:

Using the method

`zipWithIndex`

applied function`fun_index`

to each element of the list`systems`

(list of system functions) along with the corresponding index.For every pair

`(fun, index)`

a new value is calculated using the last value from`acc`

(last exact observation) for the corresponding index`index`

. Here`fun`

represents a function, and`acc.last(fun_index._2)`

gets the last value from`acc`

for index`index`

.The result of the calculation is added to

`acc`

using the operator`:+`

. This updates the list`acc`

adding a new accurate observation of the system.At the end of all iterations, a list is obtained

`positivePreciseObservations`

containing a sequence of “exact” system values *6 after each iteration in time steps`hk`

Where`k`

takes values from`-n`

before`n`

.

To generate observations for negative values `k`

perform the following steps:

For the reverse step, the formula is derived from Equation *14.

```
val negativePreciseObservations: List[List[Double]] = (1 to n).foldRight(List(yi0))((itr, acc) => {
systems.zipWithIndex.map(fun_index => acc.head(fun_index._2) - h * calkFun(fun_index._1, acc.head)) :: acc
}).init
```

Then concatenate the values for negative and positive values `k`

and perform the transposition of the list of observations.

`val preciseObservations: List[List[Double]] = (negativePreciseObservations ::: positivePreciseObservations).transpose`

## Step 3 – generating inaccurate observations.

```
val rnd = new Random()
val nonPreciseObservations: List[List[Double]] = preciseObservations.map { nums =>
nums.map(num => (1 + Math.pow(-1, rnd.nextInt()) * rnd.nextDouble() * 0.2) * num)
}
```

This step adds an error of up to 20% of the value to each value in the “exact values” list.`Math.pow(-1, rnd.nextInt())`

defines the direction of deflection (positive or negative)`rnd.nextDouble() * 0.2`

determines the amount of deviation.

## Step 4 – Implementing the Parameter Estimation Algorithm

Implementation of the parameter estimation algorithm:

```
trait SystemParameterEstimator extends StringCalculator {
def estimator(nonPreciseObservations: List[List[Double]], reverseSystems: List[String], h: Double): List[Double] = {
...
}
private def replaceVariableValues(fun: String,symbol: String , values: List[Double]): String = {
...
}
}
```

We will create a trait that will be responsible for this functionality. To work with the string calculator, do not forget to add the “extendsStringCalculator” trait.

`estimator`

– the method itself, which will return parameter estimates. He takes inaccurate observations `nonPreciseObservations`

system solution *10 `reverseSystems`

And `h`

– step.

`replaceVariableValues`

is a method that replaces `si`

to the corresponding values. In our case, this will be used to replace `yi`

And `Fi`

.

```
private def replaceVariableValues(fun: String,symbol: String , values: List[Double]): String = {
values.indices.foldLeft(fun)((acc, i) => acc.replace(s"$symbol${i + 1}", s"(${values(i)})"))
}
```

To calculate the value `n`

*9 from the available data, given that the number of steps is always equal to `2 * n + 1`

the formula will take the following form.

`val n = (nonPreciseObservations.head.size - 1 ) / 2`

Next, we will calculate the averages using *11.

`val yi0 = nonPreciseObservations.map(yi => yi.sum / yi.size)`

We will also calculate the approximate values of the functions for `yi0`

using *12.

```
val Fi0 = nonPreciseObservations.map(yi => {
(-n to n).zip(yi).map(k_y => k_y._1 * k_y._2 * h).sum / (-n to n).map(k => Math.pow(k * h, 2)).sum
})
```

Everything here is exactly according to the formula, the only thing worth noting is that the range (-n to n) is used to go through all the steps.

The last thing left is to solve the *10 system and get approximate values `bi`

.

```
reverseSystems.zipWithIndex.map( sys_index => {
val expression1: String = replaceVariableValues(sys_index._1, "F", Fi0)
val expression2: String = replaceVariableValues(expression1, "y", yi0)
calculate(expression2)
})
```

We take the equation (as a string), replace the values `Fi`

to the relevant `Fi0`

and also replace `yi`

on their average values at the zero point (k = 0).

Now the class looks like this:

```
trait SystemParameterEstimator extends StringCalculator {
def estimator(nonPreciseObservations: List[List[Double]], reverseSystems: List[String], h: Double): List[Double] = {
val n = (nonPreciseObservations.head.size - 1 ) / 2
val yi0 = nonPreciseObservations.map(yi => yi.sum / yi.size)
val Fi0 = nonPreciseObservations.map(yi => {
(-n to n).zip(yi).map(k_y => k_y._1 * k_y._2 * h).sum / (-n to n).map(k => Math.pow(k * h, 2)).sum
})
reverseSystems.zipWithIndex.map( sys_index => {
val expression1: String = replaceVariableValues(sys_index._1, "F", Fi0)
val expression2: String = replaceVariableValues(expression1, "y", yi0)
calculate(expression2)
})
}
private def replaceVariableValues(fun: String,symbol: String , values: List[Double]): String = {
values.indices.foldLeft(fun)((acc, i) => acc.replace(s"$symbol${i + 1}", s"(${values(i)})"))
}
}
```

## Last step

Let’s mix a trait `SystemParameterEstimator`

to the class where the experiment is being carried out using the keyword `with SystemParameterEstimator`

.

`object LorenzAttractor extends App with StringCalculator with ConstantFabric with BinaryOperationFabric with SystemParameterEstimator {`

Then we call the method `estimate`

on the received data and display the result.

` println(estimator(nonPreciseObservations, reverseSystems, h))`

## Result

with parameters:

```
val n = 1000
val alpha = 1.25
val h = Math.pow(n, -alpha)
val bi: List[Double] = List(2, 3, 1)
val yi0: List[Double] = List(1, 2, 1)
```

i got the result:

`List(1.9609705152167594, 3.1061318517029735, 0.9548191490247522) - b1, b2, b3 соответственно`

Well, if the results of the estimation do not differ much from the real values of the parameters, this indicates the good performance of the algorithm and the accuracy of the estimation. This may indicate that your method is effectively estimating the parameters of a system of differential equations based on inaccurate observations.

## Conclusion

I know that one example is not enough and I could do experiments on different systems with different parameters, plot frequency histograms, plot the dependence of the average deviation on the variance, and so on.

This is exactly what I’m going to do, but in the next article)

Since I have never done such work, I will ask the commentators to write on which systems it is worth doing, and what exactly is worth demonstrating in the article (frequency histograms, graphs, and so on) for a larger and more accurate study and creating a more complete picture.

## PS The project can be found on GitHub

https://github.com/PicoPicoRobotWoman/estimation_of_parameters_of_differentia_equations