Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.1
Lecture 6
Numerical Techniques : Differential
Equations - Initial Value Problems
Projectile with Air Resistance, Planetary Motion
IKC-MH.55 Scientific Computing with Python at November 17,
2023
Dr. Cem Özdo
˘
gan
Engineering Sciences Department
˙
Izmir Kâtip Çelebi University
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.2
Contents
1 Differential Equations - Initial Value Problems
Projectile Motion with Air Resistanc e
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.3
Differential Equations I
•
Most problems in t he real world are modeled with
differential equations because it is easier to see the
relationship in terms of a derivative.
•
e.g. Newton’s Law: F=Ma, d
2
s/dt
2
= a = F/M (constant
acceleration). 2
nd
order ordinary differential equation.
•
It is ordinary since it does not involve partial differentials.
•
Second order since the order of the derivative is two.
•
The solution to thi s e quation is a function,
s(t ) = (1/2)at
2
+ v
0
t + s
0
.
•
Two arbitrary constants, v
0
and s
0
, the initial values for the
velocity and posi tion.
•
The equation for s(t) al lows the computation of a nume rical
value for s, the position of the object, at any value for time,
the independent variable, t.
•
e.g. Harmonic oscillator problem in mechanics,
•
e.g. One-dimensional Schrödinger equation in quantum
mechanics,
•
e.g. One-dimensional Laplace equation in electromagnetic
theory, etc .
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.4
Differential Equations II
•
Analytical solutions of these equations are often
non-existent or ver y complicated.
•
Numerical solutions are the remedy. In terms of solution
technique, we can divide differential equations into three
groups:
1 Initial Value Problems:
In time-dependent problems, the initial state at time t=0 is given
and a solution is searched for later t values. For example, in
the following quadratic equation
d
2
y
dt
2
= f (y, y
′
, t)
two initial conditions must be given at t=0, namely y(0) and
y
′
(0) values. (N
th
order DE → N initial conditions).
2 Boundary Value Problems.
3 Eigenvalue (characteristic-value) Problems.
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.5
Projectile Motion with Air Resistance I
•
In addition to a vertical grav itational force on a 2D
projectile motion, there is also a friction force
to a certain
extent due to air res istanc e.
•
This frictional force is usually in the opposite direction to
velocity and is proportional to the square of the velocity:
~
F
r
= −kv
~
v (Drag for ce, F
D
= −(1/2)cρAv
2
~
v/|
~
v| here, c is
the drag coefficient, ρ the air density, and A the projectile’s
cross-sectional ar ea) .
If we write Newton’s 2
nd
law as
a vector in two dimensions,
m
~
a =
~
F
net
m
d
2
~
r
dt
2
= m
~
g − kv
~
v
•
and component wise (where k/m = γ):
d
2
x
dt
2
= −γ
q
v
2
x
+ v
2
y
v
x
&
dx
dt
= v
x
d
2
y
dt
2
= −g − γ
q
v
2
x
+ v
2
y
v
y
&
dy
dt
= v
y
•
Now, we have a set of equations.
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.6
Planetary Motion I
•
In the previous projectile motion example, we used the
gravitational force with the expression F = mg and
gravitational acceleration as being constant near the
Earth’s surface.
•
However, the gravitational force between masses is most
generally given by Newton’s law of universal gravitation:
F = G
m
1
m
2
r
2
Here, G = 6.6743 × 10
−11
m
3
kg
−1
s
−2
is called the
universal gravitational constant. The force is attractive and
along the direction connecting the two masses.
•
This expression should be used when studying the motion
of planets and moons.
•
Let’s study the motion of a planet (mass m )moving under
the gravitational force of the Sun (mass M). If we take the
sun at the origin, the vector expression of the force acting
on the planet would be:
~
F = −G
Mm
r
3
~
r
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.7
Planetary Motion II
•
Since the orbit of the planet will be at a plane (2D), the
position vector
~
r and accordingly the acceleration vector
~
a
would have two components as:
~
r = x
ˆ
i + y
ˆ
j
~
a =
d
2
x
dt
2
ˆ
i +
d
2
y
dt
2
ˆ
j
•
Newton’s 2
nd
law as
~
a =
~
F /m and also velocity
expressions for the x- and y-components:
d
2
x
dt
2
= −G
M
r
3
x &
dx
dt
= v
x
d
2
y
dt
2
= −G
M
r
3
y &
dy
dt
= v
y
•
Now, we have a set of equations.
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.8
Euler Method I
•
In an initial-value pro blem, the numerical solution
begins at the initial poin t and marches from there to
increasing values for the independent variable.
•
The Euler method. Describes a method that is
easy to use but is not very precise unless the step size,
the intervals for the projection of the solution, is very small.
•
Consider the following first-order differential equation:
dy
dx
= y
′
(x) = f (x, y) & y(x
0
) = y
0
(1)
•
Here x is the variable, y(x) and f(x,y) ar e real functions,
and the initial condition y
0
is a real number.
•
From the solution of this equation, we get y
1
, y
2
, . . . , y
n
values for the function at the points x
1
, x
2
, . . . , x
n
with
equal step lengths h.
•
Equations of higher order are solved by converting
them to a system of linear equations.
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.9
Euler Method II
•
The expression given by Equation 1 is written as the
forward-difference approximation at a point x
i
by Euler’s
method.
y
i+1
− y
i
h
+ O(h) = f (x
i
, y
i
)
•
If we solve this expression for y
i+1
, we get the
Euler method formula:
y
i+1
= y
i
+ hf (x
i
, y
i
) + O(h
2
)
•
This expression shows that the error in one step of Euler
method is O(h
2
). But, this error is just the local error. Over
many steps, the global error becomes O(h) (as
NO(h
2
) ≈ O(h) for N steps).
•
The method is easy to program when we know the formula
for y
′
(x)(≡ f (x
i
, y
i
)) and a starting value, y
0
= y(x
0
).
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.10
Euler Method III
•
Let’s see the application of this method on an example.
Given differential equation,
dy
dx
= x + y
•
The analytical solution of this equation is given as
y(x) = 2e
x
− x − 1. Initial condition: y(x = 0) = 1
Table: Solution of the differential equation dy /dx = x + y in the
interval [0, 1] by Eule r method.
(Example py-file: myeuler.py)
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.11
Euler Method IV
As can be seen from the table, the margin of err or is large in
the Euler method.
0.0 0.2 0.4 0.6 0.8 1.0
x
1.0
1.5
2.0
2.5
3.0
3.5
4.0
y
Approximate and Exact Solution for Simple ODE:
dy
dx
=
x
+
y
Approximate
Exact
SciPy
Figure: Solution of the differential equation dy /dx = x + y in the
interval [0, 1] by Eule r method.
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.12
Runge-Kutta Method I
•
Simple Euler method comes from using just one term from
the Taylor series for y(x) expanded about x = x
0
.
•
What if we use more terms of the Taylor series? Runge
and Kutta, developed algorithms from using more than two
ter ms of the series.
•
In the Euler method, the increment is directly from x
i
to
x
i+1
.
•
Second-order Runge-Kutta methods are obtained by using
a weighted average of two increments to y(x
0
), k
1
and k
2
.
•
Let’s take a "trial step" in the middle and then increment to
x
i+1
by using these middle x- and y-values. Two quantities
are defined here as k
1
and k
2
,
k
1
= hf (x
i
, y
i
)
k
2
= hf (x
i
+
1
2
h, y
i
+
1
2
k
1
)
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.13
Runge-Kutta Method II
•
The parameter;
•
k
1
is for the calculation at x
i
, y
i
,
•
k
2
is for a half-step away (x
i
+
1
2
h, y
i
+
1
2
k
1
) calcul ation.
•
Accordingly, the 2
nd
order Runge-Kutta formula becomes:
y
i+1
= y
i
+ k
2
+ O(h
3
)
•
In the Runge-Kutta method, the margin of error in one step
is O(h
3
) and is O(h
2
) in the entire interval.
•
It works better than the Euler method, but it comes at a
cost: f(x, y) will be calculated twice at each s tep.
•
This "trial step" technique can be taken even further.
Fourth-order Runge-Kutta (RK4) methods ar e most widely
used and are derived in similar fashion.
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.14
Runge-Kutta Method III
k
1
= f (x
i
, y
i
)
k
2
= f (x
i
+
1
2
h, y
i
+
1
2
hk
1
)
k
3
= f (x
i
+
1
2
h, y
i
+
1
2
hk
2
)
k
4
= f (x
i
+ h, y
i
+ hk
3
)
y
i+1
= y
i
+
h
6
(k
1
+ 2k
2
+ 2k
3
+ k
4
) + O(h
5
)
•
The local error term for the fourth-order Runge-Kutta
method is O(h
5
); the global error would be O(h
4
).
•
It is computationally more efficient than the (modified)
Euler method because the steps can be manyfold larger
for the same accuracy.
•
However, four evaluations of the function are required
per step rather than t wo.
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.15
Runge-Kutta Method IV
•
Let’s apply the RK4 method on the previous example.
Given differential equation,
dy
dx
= x + y
•
The analytical solution of this equation is given as
y(x) = 2e
x
− x − 1. Initial condition: y(x = 0) = 1
Table: Solution of the differential equation dy/dx = x + y in the
interval [0, 1] by 4th order Runge-Kutta method .
(Example py-file: myrungekutta.py)
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.16
Runge-Kutta Method V
As can be seen from the Table, much more sensitive results
are obtained compared to the Euler method.
0.0 0.2 0.4 0.6 0.8 1.0
x
1.0
1.5
2.0
2.5
3.0
3.5
4.0
y
Approximate and Exact Solution for Simple ODE:
dy
dx
=
x
+
y
Approximate
Exact
SciPy
Figure: Solution of the differential equation dy /dx = x + y in the
interval [0, 1] by Eule r method.
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.17
2
nd
Degree Equations & Linear Systems I
•
Any second-o rder or higher-order differential equation
can be converted into a system of first-order (linear)
equations. For example,
d
2
y
dx
2
+ A(x)
dy
dx
+ B(x)y(x) = 0
•
Let’s define two new functions for the equation, y
1
(x) and
y
2
(x):
y
1
(x) = y(x) & y
2
(x) =
dy
dx
•
With this transformation, instead of one 2
nd
order
equation, two 1
st
order equations are for med:
(1)
dy
1
dx
= y
2
(x)
(2)
dy
2
dx
= −A(x)y
2
− B(x)y
1
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.18
2
nd
Degree Equations & Linear Systems II
•
All we need to do to solve higher-order equations, even a
system of higher-order initial-value problems, is to reduce
them to a system of first-order equations.
•
Such as: One M-order equation → a system with M
first-order equations.
•
Let’s take the most general s ystem of differential equations
with M unknowns:
dy
1
dx
= f
1
(x, y
1
, . . . , y
M
) & y
1
(0) = y
10
.
.
.
.
.
. (2)
dy
M
dx
= f
M
(x, y
1
, . . . , y
M
) & y
M
(0) = y
M0
•
The next step for solving is to apply the methods (such as;
Euler, Runge-Kutta) for the 1
st
t order differential equation
to these linear system.
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.19
Projectile Motion with Air Resistance II
We had a set of equations. Two second degree and two first
degree differential equations with two unknowns.
(3)
d
2
x
dt
2
= −γ
q
v
2
x
+ v
2
y
v
x
& (1)
dx
dt
= v
x
(4)
d
2
y
dt
2
= −g − γ
q
v
2
x
+ v
2
y
v
y
& (2)
dy
dt
= v
y
•
To solve these two 2
nd
degree equations
(plus two 1
st
degree equations) given above,
we first convert them to a system of 4 1
st
degree (linear) equations.
•
To this end, let’s define the four unknowns
as follows:
•
x → y
1
•
y → y
2
•
v
x
→ y
3
•
v
y
→ y
4
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.20
Projectile Motion with Air Resistance III
•
Accordingly, the above 2
nd
degree system is written as:
(1)
dy
1
dt
= y
3
(2)
dy
2
dt
= y
4
(3)
dy
3
dt
= −γ
q
y
2
3
+ y
2
4
y
3
(4)
dy
4
dt
= −g − γ
q
y
2
3
+ y
2
4
y
4
•
When γ = 0 in this system of equations, we obtain our
usual parabolic curve y = (v
0y
/v
0x
)x − (g/2v
2
0x
)x
2
.
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.21
Projectile Motion with Air Resistance IV
To calculate the effect of air friction, let’s take the initial
conditions (t = 0) and constants (g & γ):
x
0
= y
1
(t = 0) = 0 & y
0
= y
2
(t = 0) = 0
v
0x
= y
3
(t = 0) = 6.0 & v
0y
= y
4
(t = 0) = 8.0
g = 10.0 & γ = 0.01
0 2 4 6 8 10
x
−3
−2
−1
0
1
2
3
4
y
Projectile motion with Air Resistance
Air Friction
No Air Friction
Figure: Numerical solution of projectile motion with and without air
friction. (Example py-file: airfriction.py)
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.22
Planetary Motion III
We had a set of equations. Two second degree and two first
degree differential equations with two unknowns.
(3)
d
2
x
dt
2
= −G
M
r
3
x & (1)
dx
dt
= v
x
(4)
d
2
y
dt
2
= −G
M
r
3
y & (2)
dy
dt
= v
y
•
To solve these two 2
nd
degree equations
(plus two 1
st
degree equations) given above,
we first convert them to a system of 4 1
st
degree (linear) equations.
•
To this end, let’s define the four unknowns
as follows:
•
x → y
1
•
y → y
2
•
v
x
→ y
3
•
v
y
→ y
4
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.23
Planetary Motion IV
•
Accordingly, the above 2
nd
degree system is written as:
(1)
dy
1
dt
= y
3
(2)
dy
2
dt
= y
4
(3)
dy
3
dt
= −
GM
[y
2
1
+ y
2
2
]
3/2
y
1
(3)
(4)
dy
4
dt
= −
GM
[y
2
1
+ y
2
2
]
3/2
y
2
•
For the motion of the planets, we use the astronomical unit
system. The Earth-Sun average distance would be in units
of astronomical length: 1 au ≈ 1.5 × 10
11
m. The time
taken for the Earth to go around the Sun once is 1 year (y)
as the unit of time.
•
Calculated in these units, the product of GM,
GM ≈ 40(au)
3
/y
2
Numerical Techniques:
Differential Equations -
Initial Value Problems
Dr. Cem Özdo
˘
gan
LOGIK
Differential Equations -
Initial Value Problems
Projectile Motion with Air
Resistance
Planetary Motion
Euler Method
Runge-Kutta Method
Second Degree Equations
6.24
Planetary Motion V
•
To calculate the planetary motion, let’s take the initial
conditions at time t=0 in terms of four unknowns:
x
0
= y
1
(t = 0) = 1.0 au & y
0
= y
2
(t = 0) = 0
v
0x
= y
3
(t = 0) = 0.0 & v
0y
= y
4
(t = 0) = 6.0 au/y
•
Then, also take v
0y
= y
4
(t = 0) = 8.0 au/y.
−2.5 −2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0
x
−1.0
−0.5
0.0
0.5
1.0
1.5
2.0
y
Planetary Motion
Closed
Unbounded
Figure: Numerical solution of planetary motion. There can be closed
orbits (ellipse), or solu tions going to infinity (unbounded, hyperbola)
for different velocities. (Example py-file: planetarymotion.py)
