The above classical textbook equation describes a 2nd order dynamical system.
In this formulation, `x(t)` represents
the displacement, `ω` represents the natural
frequency, and `ζ` represents the damping ratio.
To simulate this system with `IVSOLVE`, we transform it to two first order equations:

Using the named variables and initial values, we define the input to `IVSOLVE` as shown in Table 1.

A | B | C | D | |

1 | Variables with initial values | Parameters | ||
---|---|---|---|---|

2 | t | omega | 1 | |

3 | x | 1 | zeta | 0.25 |

4 | v | 0 | ||

5 | Equations | |||

6 | dx/dt | =v | ||

7 | dv/dt | =-2*zeta*omega*v-omega^2*x |

Next, we evaluate the array formula `=IVSOLVE(B6:B7, (t,x,v), {0,12})`

in array `I1:K41`
and obtain the solution shown Table 2 and plotted in Figure 1.

I | J | K | |

1 | t | x | v |

2 | 0 | 1 | 0 |

3 | 0.4 | 0.926057 | -0.35295 |

4 | 0.8 | 0.733004 | -0.59143 |

5 | 1.2 | 0.470053 | -0.70204 |

6 | 1.6 | 0.187507 | -0.69215 |

7 | 2 | -0.07064 | -0.585 |

8 | 2.4 | -0.2719 | -0.41357 |

9 | 2.8 | -0.39777 | -0.21403 |

10 | 3.2 | -0.4439 | -0.02004 |

11 | 3.6 | -0.41815 | 0.14165 |

12 | 4 | -0.33723 | 0.253773 |

13 | 4.4 | -0.22273 | 0.309249 |

14 | 4.8 | -0.09711 | 0.31042 |

15 | 5.2 | 0.019632 | 0.26696 |

16 | 5.6 | 0.112406 | 0.193178 |

17 | 6 | 0.172275 | 0.10513 |

18 | 6.4 | 0.19664 | 0.018001 |

19 | 6.8 | 0.188456 | -0.05592 |

20 | 7.2 | 0.154784 | -0.10843 |

21 | 7.6 | 0.105069 | -0.13591 |

22 | 8 | 0.049333 | -0.13896 |

23 | 8.4 | -0.00336 | -0.12157 |

24 | 8.8 | -0.04602 | -0.08994 |

25 | 9.2 | -0.07437 | -0.05117 |

26 | 9.6 | -0.08693 | -0.01211 |

27 | 10 | -0.08478 | 0.021603 |

28 | 10.4 | -0.07088 | 0.046117 |

29 | 10.8 | -0.04936 | 0.059586 |

30 | 11.2 | -0.02468 | 0.062088 |

31 | 11.6 | -0.00094 | 0.055252 |

32 | 12 | 0.018627 | 0.041749 |

The solution shows that the system is underdamped with an absolute overshoot over 0.4 at approximately a peak time of 3.3.
In this exercise, we will customize the system response by optimizing `ζ` and `ω` to achieve an absolute overshoot of exactly 0.2 at `t = 2.0`.

We define two constraints on the solution obtained in Table 2 to demand that
`x(2) = -0.2` and `v(2) = 0`, then use `NLSOLVE` to solve for `ω` and `ζ`
that staisfy these constraints. The two constraint formulas are shown in Table 3.

C | |

1 | =ARRAYVAL(J7)-(-0.2) |

2 | =ARRAYVAL(K7) |

`J7` and `K7` refere to `x(2)`
and `v(2)` values from the solution array shown in Table 2, but what is `ARRAYVAL`? `ARRAYVAL`
is one of the criterion functions that allow us to define dynamic constraints on the solution array.
In this case we are simply picking a value from the solution but in general, a criterion function can be used to specify more complex constraint on the solution.

We evaluate the array formula `=NLSOLVE(C1:C2, (w, zeta))`

in array `D1:E5` and obtain the solution
shown in Table 4 and plotted in Figure 2. The new solution has an absolute overshoot of 0.2 at `t=2` as desired.

D | E | |

1 | omega | 1.76493157 |

2 | zeta | 0.455950294 |

3 | SSERROR | 7.61331E-12 |

4 | ITRN | 11 |

5 | TIME(s) | 0.283 |

What if we wanted our peak time to be at 2.5 instead of 2. Our solution array in Table 1 does not list the output point `t=2.5`. We have two options:

We can format the output to display the solution at

`t=2.5`(See Output Format section in IVSOLVE help page).We define the constraints using another criterion function

`ODEVAL`which allows us to interpolate the solution values at any point.

Using `ODEVAL`, we modify the constraints formulas as shown in Table 5:

C | |

1 | =ODEVAL(I2:I32,x,"INTERP",2.5)-(-0.2) |

2 | =ODEVAL(I2:I32,v,"INTERP",2.5) |

The first argument to `ODEVAL`, is the time vector from the solution array in Table 2. The 2nd argument is the name of the variable
we want to operate on. The 3rd argument is the operation (here interpolation), and the 4th argument is the time value we want to interpolate at.

Using the same formula but with updated constraint formulas, `NLSOLVE` computes the new solution shown in table 6 Below
and plotted in Figure 3. Note that the peak time has shifted to 2.5.

D | E | |

1 | omega | 1.411995975 |

2 | zeta | 0.455958508 |

3 | SSERROR | 1.55938E-25 |

4 | ITRN | 8 |

5 | TIME(s) | 0.152 |