# Scilab

**Chapter 7. ****Ordinary Differential Equations**

**Chapter 7.1 Introduction**

**Scilab**is a specialized programming language for

**mathematical calculations**. Mathematical operations, e.g.

**matrix operations**, are much easier than in other languages. Therefore, he would not be himself, if he forgot about

**differential equations**. More specifically about

**ordinary differential equations**. This chapter will be less about

**Scilab**itself and more about

**mathematics**. You can skip it if you want, but I don’t recommend it.

**Scilab**will only appear in the last section. Here we will solve the differential equation

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

**recursive method**. They were also compared with the

**theoretical solution**. The method of solving ordinary differential equations in

**Scilab**works on a similar principle, only in a more sublime way. This is an

**ODE**function, that we will discuss in

**Chapter 8**.

**Chapter 7.2 Charging a capacitor as an example of the simplest differential equation**

That is, a

**first**-order differential equation, i.e. the highest derivative in this equation is the

**first**derivative

**dy/dt**, or

**y'(t)**.

**Fig. 7-1**

Charging a capacitor as an example of a first-order ordinary differential equation.

**Fig.7-1a**

**RC**system

Every electrician/electronics engineer starts with this system.

**Fig.7-1b**

Differential equation describing the process of charging a capacitor.

The

**initial state**, i.e. for

**to = 0 sec**, the capacitor voltage

**y (to)=0V**

The derivative

**y(t)**, i.e.

**dy/dt=(1/T)x(t)-(1/T)y(t)**, is the charging rate of capacitor

**C**. It is the largest for

**to=0 sec**, because capacitor

**C**is uncharged

**(y( t)=0**) and the difference

**(1/T)*x(t)-(1/T)*y(t)**is the largest. Then the

**capacitor**will charge a little and this difference will decrease. That is, the charging speed, otherwise the

**dy/dt**slope will also decrease. Until we reach the steady state when

**dy/dt=0**. Otherwise,

**y(t)**will become a horizontal line.

**Chapter 7.3 Generalization of a first-order differential equation.**

**Fig.7-2**A more generalized differential equation

**Fig. 7-2a**

Generalization of the equation

**Fig.7-1a**

This is an example of a

**first-order linear differential equation**with the initial state

**to=0, y(to)=0**.

Parameters

**a**and

**b**are arbitrary constants, e.g.

**a=1 b=1**

**Note:**

I am an automation engineer and that is why I associate

**differential equations**most often with

**time**functions

**y(t)**, and less often with the

**shift**function

**y(x)**. For a mathematician it has no meaning.

**Fig. 7-2b**

Diagram with an integrating system for solving the equation in

**Fig. 7-2a**.

Back in the

**1950s…1960s**, computers were called

**analog**and

**digital**mathematical machines.

**Digital machines**are simply

**computers**, and

**analog**machines are systems with

**integrating elements**, somewhat forgotten nowadays, as above. In the

**red rectangle**there is a diagram performing the mathematical operation

**ax(t)-by(t)**, or generally

**f(t,y)**. After the integrating circuit there is voltage

**y(t)**. In the initial state,

**to=0**often has the value

**y(to)=0**. As for any integrating system, “after” is

**y(t)**, and “before” is the derivative

**dy/dt**. The system will begin solving the differential equation starting from

**to = 0**and

**y(to) = 0**, producing a waveform as shown in Figure

**7-2c**.

**Fig. 7-2c**

Solving the differential equation in

**Fig. 7-2a**by the system in

**Fig. 7-2b**.

**a=1, b=1 y(to)=0**and when

**x(t)=1(t)**is a

**unit step signal**.

The analytical solution, or “pencil by a mathematician”, is

**y(t)=1-exp(-t)**. If you don’t believe it, substitute

**y(t)=1-exp(-t)**into the differential equation

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

**Chapter 7.4 The most general first-order differential equation**

**Fig.7-3a**

The most general first-order differential equation

On the right side of

**dy/dt**is an any expression

**f(t,y)**, where

**t**is time and

**y**is a function of time

**y(t)**.

In simple cases it may be e.g.

**f(t,y)=y(t)**.

**Fig.7-3b**

Diagram with an integrating system for solving the equation in

**Fig. 7-3a**.

**Chapter 7.5 Examples of first-order differential equations**

The expression

**f(t,y)**in the equation

**Fig. 7-3a**is not very telling. Some specific examples will be useful

.

**Fig.7-4a**

First order linear differential equation

The signal **x(t)=1(t)** is a **unit step**, i.e. **1(t)=0** **for t<0** and** 1(t)=1 for t>0**.

Here and in the remaining drawings I did not show the initial state, e.g. **to=0** and **y(to)=0**.**Fig.7-4b**

The differential equation is **non-linear** because it contains **2y²(t)**.

The function** t(t)** is the so-called **“ramp”**—>**t(t)=0 for t<0 and t(t)=t for t>0**.**Fig.7-4c**

The differential equation is **non-linear**, because **y(t)** is “under the square root”. Signal **1(t-3)** is a step delayed by **3** seconds.**Fig.7-4d**

Another example of a nonlinear differential equation.

I will also add, that the input signal **x(t)** may be non-linear, e.g. **x(t)=t²(t)** or steps **1(t)**, **1(t-3)**, although e.g. **dy/dt=3y(t)+t²(t)** is a linear differential equation!**Chapter 7.6 How the model in Fig. 7-3b solves a differential equation?**

We will analyze the system solving the differential equation in **Fig. 7-3b**. There, time worked continuosly. We will quantize it, i.e. we will examine its states in time up to** t4** with a time interval of **Δt**. We know that the approximate value of the **dy/dt** derivative is the so-called difference quotient **Δy/Δt**. That is, **dy/dt≈Δy/Δt**.

It turns out that knowing the initial state **y(to)** we can calculate**Δ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**.

In this way, we will calculate all values of

**y(ti)**from the initial time

**to=0**to the final time here

**tk=10sec**. So, knowing the initial state

**y(to)=yo**, we will solve the

**first**-order

**differential equation dy/dt=f(t,y)**. This is, of course, an approximate solution, the smaller

**Δt**is, the more accurate it is.

**Fig.7-5**

The state of the model in **Fig.7-3a** at time from **t to t4****Fig. 7-5a**

The output of the **integrating** element is the initial state **y(to)=y0**. The **“pseudo-derivative”,** or **difference quotient Δyo/Δto=f(t0,y0)** will be calculated immediately as a certain constant value. During the time **Δto** from **to** **t1**, the **integrating** unit “integrates” with the speed **Δyo/Δto** from the value **y0** to y**0+Δyo**. **Fig. 7-5b**

This is a photograph of **Fig.7-5a** at its final moment **t1**. So **y1=y0+Δy0**. This is also the initial state in **Fig. 7-5c**.**Fig. 7-5c**

The initial moment **t1** is the same as in **Fig. 7-5b**. This is actually the same state, except that at time** t1** a new difference quotient **Δy1/Δt1=f(t1,y1)** is calculated as a new **“pseudo-derivative”**. During the time **Δt1** from **t1** to **t2**, the **integrating** unit “integrates” at the rate **Δy1/Δt1** from the value **y1** to **y1+Δy1**. **Fig. 7-5d**

This is a photograph of **Fig.7-5c** at its final moment **t2**. So **y2=y1+Δy1**. This is also the initial state in **Fig. 7-5e**.**Fig. 7-5e**

Description similar to **Fig. 7-5c**.**Fig. 7-5f**

Description similar to **Fig. 7-5d**.**Fig. 7-5g**

Description similar to **Fig. 7-5e**.**Fig. 7-5f**

Description similar to **Fig. 7-5d**.

The next states for **t5, t6, …tk** are also calculated in a similar, i.e. recursive, way.**Chapter 7.7 How to solve a specific equation, e.g. dy/dt=y(t)-1(t)**

That is, the solution for the equation in **Fig. 7-4a** as a special case in **Fig. 7-2a**, for **a=1, b=1 and x(t)=1(t)**. The solution to this equation is **Fig. 7-2c**. Will the approximate solution based on **Chapter 7.6** confirm this?

**Fig.7-6**

An integrator as a machine for solving the equation **dy/dt=1(t)-y(t)****Fig.7-6a**

Differential equation with zero initial conditions**Fig.7-6b**

Integrator as a machine solving the equation **dy/dt=1(t)-y(t)****Fig.7-6c****y(t)=1-exp(-t)** as an exact solution (“in pencil”) of the equation **dy/dt=1(t)-y(t)**

I suggest substituting **y(t)** into the differential equation and checking. Remember that derivative** [exp(-t)]’=-exp(-t)**

Now we will solve this equation numerically, assuming **dy/dt≈ Δy/Δt**. The smaller the time increment **Δt,** the more precise it is.

**Fig.7-7**

Numerical solutions of the differential equation **y'(t)=1(t)-y(t)****Fig.7-7a**

Exact differential equation–>** y'(t)=1(t)-y(t)**

Difference equation as approximate differential equation–>**Δy/Δt=1(t)-y(t)****Fig.7-7b**

How was the recursive equation **y(n+1)=y(n)(1-Δt)+Δt** obtained from the difference equation?

My grandson took his **8th grade exam** yesterday. If the differential quotient were treated as an ordinary quotient, having nothing to do with differential calculus, it would arrive at the equation in the red box. Note that for a given increment of** Δt,** the coefficient **(1-Δt)** is constant, which makes the calculations very easy, e.g. for **Δt=0.1 (1-Δt)=0.9**.**Fig. 7-7c**

The values of **y(2),y(3),y(4)…y(n+1)** calculated based on the equation **Fig. 7-7b** and the known values of **y(1)** and **Δt**. **y(1)=0** as the known initial value of the equation. **y(2)=Δt** always**y(3)=y(2)*(1-Δt)+Δt** is where the recursion begins.

Note that for a given **Δt**, e.g. **Δt=0.1sec**–>**1-0.1=0.9**, this recursion is very simple **y(n)=y(n-1)*0.9+0.9****Fig. 7-7d**

Next **y(n)** when **Δt=1sec**

Check for several **n.** Not very carefully, because **Δt=1sec** is large compared to the range **t=10sec**! If we increased e.g. **Δt=2sec** it would be even worse! Check. It goes with intuition. The smaller **Δt**, the more similar the difference quotient **Δy/Δt** is to the derivative **dy/dt**!**Fig. 7-7e**

Next **y(n)** when **Δt=0.5sec****Fig. 7-7f**

Next** y(n)** when **Δt=0.25sec****Fig. 7-7g**

Next **y(n)** when **Δt=0.1 sec****Conclusion**

We expect that by decreasing **Δt** the response **y(n)** approaches the ideal solution **Fig. 7-6c**.**Chapter 7.8 Solving the above equation using Scilab.****Chapter 7.8.1 Analysis of the first part of the program**

For now, we will analyze the first part of the program in which we will create and check the increase **Δt** of the difference quotient as well as the **t** vectors and the **step unit**.**7-1 program**

clc;clear;clf jtmax=10 Experiment_time=10 Δt=Experiment_time/jtmax disp('Δt=',Δt) t=0:1:jtmax disp('t=',t) step=t for j=1:jtmax+1 step(j)=1 end disp('step=',step)

Let me remind you what each instruction does**clc;clear;clf**

Standard cleaning of the console screen, variables and graphs**jtmax=10**

The expected number of calculations here is **10**, but in fact it will be** jtmax + 1**, i.e. **11**.

Usually **jtmax** is larger, e.g. **100, 1000…****Experiment_time=10**

Our differential equation covers time **t=10 sec**.

Otherwise, the domain of the differential equation is time** t = 10 sec**.**Δt=Experiment_time/jtmax****Δt** time increment calculated for the **difference quotient**.**disp(‘Δt=’,Δt)** **Δt** display **t=0:1:jtmax**

Vector **t** with discrete times **t=[0,1…jtmax]** on which the vectors will be defined **disp(‘t=’,t)**

Display vector **t** **step=t**

Definition of a new **step** vector. For now, the **step** is the same as **t****for j=1:jtmax+1****step(j)=1****end**

Conversion of the **step** vector from** “time”** to** “step”** (with “ones”). Here you may wonder why it is there**for j=1:jtmax+1** and not **for j=0:jtmax**? Because the vectors** t** and **step** do not start from** t(0)**, **step(0)**, but from **t(1)**, **step(1)**. Because that’s the typical of all vectors in **Scilab**.

disp(‘**step**=’,**step**)**step** vector display

In the animation **Fig. 4-2** in **chapter 4** I showed how to execute the program in **Scilab** and there should be no problems with executing the above.

Follow this program. You will get something like this.**“Δt=”**

**1.**// Since the experiment time

**t=10 sec**and

**jtmax=10**then

**Δt=1sec**is quite large. Usually

**Δt**is much smaller.

**“t=”**

**0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.**

**“step=”**

**1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.**

Everything, i.e.

**Δt**and v

**ectors**, are as we predicted.

**Chapter 7.8.2 Analysis of the entire program**

Now that we have calculated

**jtmax**based on it

**-Δt**– time increment

**-t**time vector with time units from

**0**to

**10 sec**

**-step**vector with

**jtmax**time units

then we can compare the differential solution

**y(j)**from

**Fig. 7-7d…7-7g**with the ideal solution from

**Fig. 7-6c**on the graph. We expect that as

**Δt**decreases, the differential solution will approach the

**ideal**one.

**Program 7-2**

clc;clear;clf jtmax=10 Experiment_time=10 Δt=Experiment_time/jtmax t=0:1:jtmax step=t for j=1:jtmax+1 step(j)=1 end//up to this point, description inchapter 7.8.1ideal=step for j=1:jtmax+1//see first explanation ideal(j)=(1-exp((1-j)*Δt)) end real=step real(1)=0 for j=2:jtmax+1 //see second and third explanation real(j)=Δt+(1-Δt)*real(j-1) end plot(t,step,'o') plot(t,step) plot(t,ideal,'o') plot(t,ideal) plot(t,real,'o') plot(t,real,'r')

The **first explanation** is when, for example, **jtmax=10–>Δt=1sec**

That is, we calculate the ideal graph of the solution of the equation in **Fig. 7-6c** for **Δt = 1 sec****j=1**–>ideal(1)=y(1)=1-exp(1-1)*1=**0****j=2**–>ideal(2)=y(2)=1-exp(1-2)*1=**0.6321…****j=3**–>ideal(3)=y(3)=1-exp(1-3)*1=**0.8646…****j=4**–>ideal(4)=y(4)=1-exp(1-4)*1=**0.9502…**

…**j=10**–>ideal(10)=y(10)=1-exp(1-10)*1=**0.9998…**

We see that the calculations of the ideal solution coincide with **Fig.7-1c**

**The second explanation** when **jtmax=10–>Δt=1sec**

That is, we calculate the course of the approximate solution of the equation in **Fig. 7-7d** for** Δt = 1 sec****j=1**–>real(1)=y(1)**=0****j=2**–>real(2)=y(2)=Δt=**1****j=3**–>real(3)=y(3)=y(2)*(1-1)+Δt=Δt=**1** **j=4**–>real(4)=y(4)=y(3)*(1-1)+Δt=Δt=**1**

…**j=11**–>real(11)=y(11)=y(10)*(1-1)+Δt=Δt=**1**

**Fig.7-8**

Graph of the

**ideal**and

**iterative**solution for

**Δt=1sec**

The step unit

**1(t)**perfectly coincides with the “top” of the drawing. The iterative approximation is the

**upper dots**and the

**lower ones**are the perfect solution. You can see the difference between the

**ideal**and

**iterative**solution better than in the

**“numbers”**in the

**first**and

**second**explanations. One more thing. In the graph, the circles were connected by straight lines, making the chart look like a

**real**one.

Strange, the approximate course is almost

**rectangular**. But this is the

**first**and very

**rough approximation**when

**Δt=1 sec**. But if I hired a good lawyer, he would prove that an

**exponential**chart is almost a

**step**!

It will definitely be better when

**Δt=0.1**. Let’s check

**Third** explanation when** jtmax=100–>Δt=0.1 sec**

That is, we calculate the chart of the actual **solution** of the equation in **Fig. 7-7g** for **Δt = 0.1 sec****j=1**–>real(1)=y(1)=**0** **j=2**–>real(2)=y(2)=Δt=0.1–>**=0.9****j=3**–>real(3)=y(3)=y(2)*(1-01)=y(2)*0.9=**0.1 9****j=4**–>real(4)=y(4)=y(3)*0.9=**0.271****j=5**—>real(5)=y(5)=y(4)*0.9+Δt=Δt=**0.3419**

…**j=101**–>real(101)=y(101)=y(100)*0.9+Δt=Δt=**0.9999734…Fig.7-9**Graph of the

**ideal**and

**iterative**solution for

**Δt=0.1 sec**

The chart is clearer because I gave up the dots. There would be

**101**of them and it would make the chart less transparent. It is evident that the

**iterative**graph came close to

**perfect**. Let’s check this for

**jtmax = 250**, when

**Δt = 0.04 sec**. It will definitely be better, but by how much? We will give up the computational part and only provide a

**graph**.

**Fig.7-10**Graph of the

**ideal**and

**iterative solution**for

**Δt=0.04 sec**

As we expected, the

**ideal**graph almost coincides with the

**iterative**one. It also coincides with intuition. The

**smaller Δt**, the

**greater**the accuracy of solving the

**differential equation.**Will it be like this until the end? There are limits to reducing

**Δt**.

**First**, the computer itself calculates with

**finite**accuracy.

**Secondly**. The smaller

**Δt**, the larger

**jtmax**, the more

**iterations**. And during

**iteration**, errors accumulate. Maybe this will end our considerations about the accuracy of mathematical modeling. But we came to one basic conclusion. There is some optimal

**Δt**. Not too big, not too small.