# Scilab

**Chapter 8. ****Function for solving Ordinary Differential Equations**

**Chapter 8.1 Introduction**

**Fig.8-1**

How does the integrator solve the first-order differential equation

**y'(t)=f(t,y)?**

The system with an **integrating unit** solves the differential equation **y'(t)=f(t,y)** when initial conditions **y(to)=yo**. In other words, the entire differential equation is “enclosed” in the **red** box in **Fig.8-1a**. The comment will also apply to further diagrams with integrating circuits. Apparently it follows from the system itself, that at every moment** t** there is clearly **y'(t)=f(t,y)**. Especially since the output of each integrating system is **y(t)** and the input is the derivative **y'(t)**. But where is the beginning and end of events? A snake eating its own tail! Just some doubts. They should disappear when analyzing **Fig.8-1b**. At least here you can see what is the beginning. Namely, the initial state **y(t0)=yo**. So **for to=0sec**, there is a voltage **yo** at the output of the integrator. Most often, **yo=0**. **Fig.8-1b** shows, that for the initial state **(to,yo)** the difference quotient **Δy0/Δto≈f(t0,yo)** can be calculated, i.e.** y1=y0+Δy0**. And what’s next? We will again calculate the new **Δy1/Δt1≈f(t1,y1)**, i.e. **y2=y1+Δy1**. Etc. This is how we will get the sequence of the output waveform **y0,y1,y2…yn**. This is an approximate solution to the differential equation **dy/dt=f(t,y)** with initial conditions **y(to)=yo**.**Note:**

The analog circuit from the **1950s** is the ancestor of the **ODE** function – a basic tool for solving **first-order ordinary differential equations**. And later, as it turns out, also of **any** order.

**Section 8.2 Solving the y'(t)=1-y(t) equation using ODE****ODE** – a function that solves arbitrary **ordinary differential equation of first order****y'(t)=f(t,y)**

at the initial conditions**y(to)=y(o)**

This is a numerical implementation of the integrating system from **Fig. 8-1** for any function **f(t,y)** with initial conditions **y(to)=f(to).**

The numerical solution for the specific function **f(t,y)=1(t)-y(t)** is given in **Chapter. 7** in **Fig.7-7** and this is implemented by program **7-2**.

In a similar way, only for any **f(t,y)**, the **ODE** function works.

Let’s check the operation of **ODE** for **f(t,y)=1(t)-y(t)** in the following program**P8-1** program with **ODE** function without comments

clc; clear; clf(); function ydash=f(t, y) ydash=1-y endfunction t=[0:0.1:10] t0=0 y0=0 y=ode(y0,t0,t,f) plot(t,y,'--or')

The same program with comments. Analyze it

clc;//Console cleaning clear;//Data cleaning clf();//Graph cleaning function ydash=f(t, y)//Here startsODEfunction as general ydash=f(t,y) ydash=1-y//Here is specific differential equation to solve endfunction//Here endsODEfunction t=[0:0.1:10]//Time range y0=0//initial stateyowhent=0y=ode(y0,t0,t,f)//ODEzalculates the solution based onto,yoplot(t,y,'--or')//ODEcreates graphy=f(t)

**Note:****y’** is also known as **ydash–>dash**, an apostrophe after–>**y** a derivative symbol.

Execute program **P8-1**. The **animation** below will be helpful.

**Fig.8-2**How to execute program

**P8-1**?

You see the final animation result with the calculated

**y(t)**graph, which is the solution to the differential equation

**y'(t)=1-y(t)**

**1.**Place the Scilab

**Course**and Scilab

**Console**windows next to each other.

**2.**Call up the

**Scinotes**editor (you have skipped the last edit file here)

**3.**Copy the

**P8-1**program from the Scilab

**course**

**4.**Type

**P8-1**into the Scinotes

**editor**

**5.**Execute wthout echo. Here

**Scilab**tells you to save it somewhere.

**6**. Save the

**P8-1**program anywhere, under any name.

**7.**A graph of

**y(t)**has been made as a solution to the equation

**y'(t)=1-y(t)**with initial conditions

**y(0)=0**.

The result is the graph in

**Fig. 8-2**whose

**“circles”**correspond to the intervals

**Δ=0.1sec**from the instruction

**t=[0:0.1:1]**.

I suggest you run the

**P8-1**program yourself. Just change the instruction

**plot(t,y,’–or’)**to

**plot(t,y,’r’)**. Then the graph will be continuous.

**Fig. 8-3**

Solving the differential equation

**y'(t)=1-y(t)**as a graph of

**y(t)**.

In the next chapter you will learn that

**y(t)**is the response of the i

**nertial**unit

**K=1**and

**T=1 sec**to a

**unit step**.

**Chapter 8.3 Solving the differential equation y'(t)=1.5-2y*(t) using ODE**

The flagship of the portal is **Automatics**, in which all experiments were performed in **Scilab**, or more precisely in its **XCOS** application. Another thing is that the reader is not dealing with **Scilab** here. And who has it? The **author** of the **course** who, for the convenience of readers, recorded animations of exercises performed with **Scilab** using **ActivePresenter **applicaton. Let’s go back to the **Automatics **course. At the beginning of the course you learn the so-called dynamic elements. **Proportional, Inertial, Oscillatory.**.. etc. Their dynamics is described by the so-called **transmittance** **G(s)**. The **inertial** unit is shown below.**Fig. 8-4**

Testing the **Inertial Unit** with the **ODE** function**Fig. 8-4a**

Inertial Unit**Fig. 8-4b**

Transmittance **G(s)** of the** Inertial Unit****K** gain**T** time constant**Fig. 8-4c**

Differential equation describing the transmittance of the **Inertial Unit**.**Note:**

The relationship between the transmittance **G(s)** and the **differential equation** is shown in the course Automatics **Chapter 15.3****Fig. 8-4d**

Differential equation of the **Inertial Unit** for **K=0.75** and **T=0.5sec**. When the input signal **x(t)** is a unit step **1(t)**, we will obtain the equation from the chapter title, i.e.** y'(t)=1.5-2y*(t)**.

Let’s examine this unit with **Scilab** using the **ODE** function. The interpretation from the **graph** of the **gain** parameters **K** and the time constant** T** is obvious. **P8-2** program with **ODE** function

clc; clear; clf(); function ydash=f1(t, y) ydash=1.5-2*y endfunction t=[0:0.1:10] t0=0 y0=0 y=ode(y0,t0,t,f1) plot(t,y,) t0=0 y0=0 y=1 plot(t,y,)

Description of the **P8-2** program with the **ODE** function

How is it different from **P8-1**?**1.** instruction **ydash=1.5-2*y** instead of **ydash=1-y****2.** **plot(t,y,)** instead of **plot(t,y,’–or’)** –> solid blue graph.**3.** An additional sequence of instructions that plots the function **y=1**. This ensures that the graphs in **Fig. 8-5** and **Fig. 8-3** have the same scales.

t0=0 y0=0 y=1 plot(t,y,)

Execute program **P8-3.** You will get the graph in **Fig. 8-5**

**Fig. 8-5**

Solving the differential equation

**y'(t)=1.5-2*y(t)**as a graph of

**y(t)**.

Referring to

**Fig. 8-4**, this is the response of the

**Inertial Unit**with parameters

**K = 0.75**and

**T = 0.5 sec**per step unit.

The interpretation of the

**K**and

**T**parameters is obvious. Look again at

**Fig. 8-3**with the

**Inertial Unit**

**K=1**and

**T=1sec**.

Chapter 8.4 Solving the differential equation y”(t)=0.5-0.125*y'(t)-0.25*y(t) using ODE

Chapter 8.4 Solving the differential equation y”(t)=0.5-0.125*y'(t)-0.25*y(t) using ODE

**Section 8.4.1 General information about the above equation and its relationship to the Oscillation Unit.**

As it turns out, we are dealing here with a response to a

**step unit**of an

**Oscillation Unit**with parameters

**K=2, T=2 sec**and

**q=0.5**.

**Fig.8-6**

The **Oscillation Unit** and its differential equation**Fig.8-6a**

Transmittance **G(s)** of this unit

At first glance it is not obvious that **G(s)** is an **Oscillatory Unit**. Or maybe it’s **2-Inertial Unit**?**Fig.8-6b**

The differential equation directly related to **G(s)** in **Fig. 8-6a**.

Even if you don’t know what transmittance is, try to find the dependence on **G(s)**.**Fig.8-6c**.

The differential equation when **x(t)** is **step unit** i.e. **x(t)=1(t)**. So **0.5x(t)=0.5*1(t)=0.5**.

Because **1(t)=0** for** t<0** and **1** for **t>1**, and the range of the equation is for **t>0**.**Fig.8-6d**

Exactly the same transmittance as **Fig.8-6a**, only written differently.

But from it you can obtain the parameters of the oscillating element **K = 2**, **T = 2 sec** and **q = 0.25**.

See the** Automatics Course**, **chapter 7** **Fig. 7-3**.**Fig.8-6e**

A system with **2** **Integrating Units** describing the equation in **Fig. 8-6c****Note** that it solves the equation **y”(t)=0.5-0.125y(t)-0.25y(t)**, similarly to the system with **one** **Integrating Unit **from **Fig. 7-6 chapter. 7** solves the equation **y'(t)=1(t)-y(t)**. The **second-degree differential equation** is also included in the **red** box, as in **Fig.8-1a**.

**Chapter 8.4.2 How does ODE solve the equation y”(t)=0.5-0.125*y(t)-0.25*y(t)?**

I will not analyze in detail how to approximately solve the **2nd degree ordinary differential** equation** y”=x(t)-a*y'(t)-b*y(t)**. Or more generally **y”(t)=f[t ,y'(t),y(t)]** .

For the** first degree** equation** y'(t)=f[t,y(t)]** we replaced it with the approximate difference equation **Δy/Δt=f[t,y(t)]** with initial conditions **y(to)=yo** and then:**Δyo/Δt=f(to,yo) –>Δyo=f(to,yo)Δt**

So we will calculate **y1=y0+Δyo=y0+f(to,yo)Δt**

Similarly, knowing** y1**, we will calculate **y2** (because we will calculate **Δy1/Δt1=f(t1,y1)** etc…

That is**y2=y1+f(t1,y1)Δt****y3=y2+f(t2,y2)Δt****y4=y3+f(t3,y3)Δt**

We assume that the time increments **Δt** are equal to –>**Δt0=Δt2=Δt3=Δt4=…=Δt**.

Note that based on the previous state **y(n)** we will calculate the next state** y(n+1)**. This is the so-called recursive equation We know the first, initial state **y(0)=yo**.

The same will happen with the **2nd** degree equation. But here we have **2** initial states, the function **y(to)** and its derivative **y'(to)**. Here we will also replace the **first** derivative with the difference quotient **Δyo/Δt=f(to,yo)** and the **second** derivative (“derivative from derivative”) with the difference quotient (“increment from increment”) **Δ(Δy/Δt)/Δt**. Without going into details, this will lead us to a “double” recursive equation **y(n+2)=f[y(n+1)]** and **y(n+1)=f[y(n)]**. So, having **2** initial conditions **y(to)** and** y'(to)** we will calculate all **y(2),y(3),y(4)…y(n),y(n+1)**. Wait a minute. Where’s **y(0),y(1)**? This can be easily calculated from the initial values **y(to)=y(o)** and **y'(to)=y'(o)**. Conclusion, we will approximately solve a **2nd** order **differential equation**.

Dear reader, maybe this is too general and you didn’t understand everything. Do not worry. This is a task for the **ODE** function. You just need to enter the appropriate parameters. But the most important thing is, that we will replace **one **second-order equation with **two **first-order equations.

**Chapter 8.4.3 How does ODE treat the equation y”(t)=0.5-0.125*y(t)-0.25*y(t)?**

First of all, it turns **one** second-degree differential equation into **two** first-degree equations. How? And as below.**Fig.8-7**The

**oscillation**unit and its

**differential**equation as

**two**first-degree

**equations**

**Fig.8-7a**

**Oscillation**Unit

So

**Fig.8-6d**is equivalent to

**Fig.8-6a**

**Fig.8-7b**

Differential equation describing the

**Oscillation**Unit

**Fig.8-7c**

As above but with

**two**equations

**Fig.8-7d**

A version of

**two**first-degree equations instead of

**one**second-degree equation from

**Fig.8-7b**

**Fig.8-7e**

An analog circuit with

**two**

**integrators**solving the differential circuit in

**Fig. 8-7d**

I will add that it will give the same answer as the system in Fig. 8-6e, which is probably obvious.

**Conclusion:**

**One**

**second**-degree differential equation can be converted into a system of

**two**

**first**-degree differential equations

**General conclusion**:

**One**

**third**-degree differential equation can be converted into a system of

**three first**-degree differential equations.

**e.t.c**…

And what about

**ODE**?

**P8-3**program

clc; clear; clf(); function ydash=f(t, y) ydash=[y(2);0.5-0.125*y(2)-0.25*y(1)]//after the symbol ";" we have entered the right side of the equation y2'(t)=...from Fig.8-7d endfunction t=[0:0.1:100]//range 0<t=<=100 with intervals Δt=0.1 t0=0//initial values t0=0sek y0=[0;0]///initial values i.e.y2(0)=0andy1(0)=0y=ode(y0,t0,t,f)//functionODE callplot(t,y(1,:)) xgrid(1,1,3)

Program **P8-3** description **ydash=[y(2);0.5-0.125y(2)-0.25*y(1)]**

This is a record of** one** second-degree differential equation with one unknown **y(t)** from **Fig.8-7b** in the version of **two** first-degree equations from **Fig.8-7d**.**Note:**

In program **P8-1** it was **ydash=1-y**. Are there any square bracket symbols here? Because they concern two variables **y2(t)** and** y1(t)** from **Fig. 8-7c**. Written as y**1′(t)** and **y2′(t)** from **Fig. 8-7d**.

That is, the vector **[y1′(t),y2′(t)]** or the vector **[y(2);0.5-0.125y(2)-0.25*y(1)]**. Here we marked **y2(t)** with the symbol **y(2)** and **y1(t)** with the symbol **y(1)**. Do not associate it with any recursion!!!

Copy and execute this program in the **Scinotes** application as in the animation **Fig.8-2**.

You will get a chart like this.** **

**Fig. 8-8**

Graph

**y(t)**as a solution to the equation

**Fig.8-7b**with initial conditions

**y(0)=0**and

**y'(0)**

There may be some doubts about the derivative in

**Fig. 8-8**.

**y'(0)≠0**–>is not zero! ! Yes, but it should be

**y'(0)=0**with the note – for a homogeneous equation, i.e. for

**y'(t)=-0.125y'(t)-0.25y(t)**. And not the one with

**0.5**as in

**Fig.8-7b**.

**Note:**

You can read more about this

**oscillation unit**in the course

**“Automatics”**

**Chapter 7**

**Chapter 8.5 Solving the nonlinear differential equation y”(t)=0.5-0.125*y'(t)-0.25*√y(t) using ODE****Fig.8-9**

Nonlinear **Oscillation Unit**

Why non-linear? Because instead of **0.25*y** (as in **Fig. 8-7b**), in **Fig. 8-9a** there is **0.25*√y** (0.25*square root of y”)**Fig.8-9a**

Due to non-linearity, the dynamic unit cannot be described as the transmittance **G(s)** as in **Fig. 8-7a**. But its dynamics more generally as a differential equation.**Fig.8-9b**

The substitution is analogous to **Fig. 8-7c**. Nonlinearity doesn’t matter.**Fig.8-9c**

Similar to **Fig.8-7d****Fig.8-9d**

Similar to **Fig.8-7e**. The difference is only in the **root unit**. The whole thing is solved by the differential equation from the **title** of the **chapter**.

The digital version of **Fig.8-9d** is the **P8-4** program below with the **ODE** function. Comments are probably unnecessary here. The only difference is the root function **sqrt(y(1))** used in **ODE**

**Program P8-4**

clc; clear; clf(); function ydash=f(t, y) ydash=[y(2);0.5-0.125*y(2)-0.25*sqrt(y(1))] endfunction t=[0:0.1:100] t0=0 y0=[0;0] y=ode(y0,t0,t,f) plot(t,y(1,:)) xgrid(1,1,3)

Copy and execute it similarly to program** P8-4. **

**Fig. 8-10**

Solution of a nonlinear differential equation

Compared to

**Fig.8-8**, the response is more muted and lazy. This is the effect of signal square rooting.