# First Order Models

Let's start by creating a simple first order model. The equations for our model will be:

\frac{dx}{dt} = -x

Note, we will do this without using the "component oriented" features of ModelingToolkit just so you can see how a very basic model is constructed.

Let's start by using the following packages:

using ModelingToolkit
using DifferentialEquations
using Plots

These using statements make the contents of these packages (at least the bits that are explicitly exported) available for us to use in our file.

Now, we need to create a few definitions up front. For example, we need to declare the variables in our system, e.g.,

@variables t x(t)

Anything that starts with @ in Julia is a macro. For now, don't worry about what that means precisely. The important point is that this is how you let ModelingToolkit know that both t and x are variables and that x is a function of t.

You might wonder why it is necessary to state that x is a function of t. Part of the reason is that MTK supports PDEs as well as ODEs and this is a way of telling MTK that x is only a function of t and not a function of any other spatial variables.

In order to write our differential equation, we need to introduce an operator that represents the derivative of something. We do this as follows:

D = Differential(t)

This means that if we apply D to a variable, we should treat the result as the derivative of that variable. We will use this to represent derivatives in our equations.

Speaking of equations, we can now create an ODESystem from a set of equations:

@named sys = ODESystem([D(x) ~ -x], t)

Don't worry so much right now about what the @named macro does, that will disappear soon anyway. The main thing to notice here is that we create an ODESystem by passing in a vector of equations. You should read D(x) ~ -x as "The derivative of x with respect to time equals -x". Note that equations use a ~ to represent "equality", not =. Not that D doesn't actually take the derivative of x here, how could it? We don't (yet) know how x varies with t. Instead, D(x) is used to symbolically represent the derivative of x with respect to t.

We can now take this system of equations, represented by the variable sys, and we can formulate a complete initial value problem (an ODEProblem) as follows:

prob = ODEProblem(structural_simplify(sys), [x => 10.0], (0, 10), [], jac = true)

...where the first argument, structural_simplify(sys), is our system of equations (after they've been simplified), the second argument is the set of initial conditions, the third argument is that interval of time we wish to solve these equations over, the fourth argument contains parameter values (empty here, since we don't have parameters) and finally we have settings specified by keyword arguments.

In this case, our problem is to solve the equation D(x) ~ -x over the interval from 0 seconds to 10 seconds, assuming the initial value for x is 10.0.

Now that we have a problem, we can solve the problem by running:

sol = solve(prob)

Once we have a solution in hand, we can plot the result with:

display(plot(sol, idxs=[x]))

We can combine all of this together into a single file, i.e.,

using ModelingToolkit
using DifferentialEquations
using Plots

@variables t x(t)

D = Differential(t)

@named sys = ODESystem([D(x) ~ -x], t)

prob = ODEProblem(structural_simplify(sys), [x => 10.0], (0, 10), [], jac = true)

sol = solve(prob)
display(plot(sol, idxs=[x]))

This single file contains all the package dependencies, declarations of variables, equations systems, simulation problem setup, solution and plotting all in a single script.

If you run this file from the command line, e.g., runnings

julia first_order.jl

...or use the "Julia: Run File in New Process" command in Visual Studio Code, you won't see any output. In order to render the plot, this script can be run in Visual Studio Code by selecting the "Julia: Execute File in REPL" command or you can add the -i command flag, e.g.,

julia -i first_order.jl

Note, it will take some time for the plot to appear. This is because Julia is performing a "just in time" compilation of the model. This pre-compilation only needs to be done once per sessions. So, for example, if you've already run that script above you can do the following very quickly:

prob = ODEProblem(structural_simplify(sys), [x => 10.0], (0, 100), [], jac = true)
sol = solve(prob)

display(plot(sol, idxs=[x]))

...either by typing this into the Julia REPL or by re-evaluating these lines in the VS Code (or re-evaluating the entire file).

The point is that once this model is compiled, simple changes (like the end time) will be compiled and simulated very quickly.

The results for this particular model should look like this: