Math 157: Intro to Numerical Solutions to ODE

Yiwei Sang, Winter 2018

Ordinary Differential Equation (ODE)

• An ordinary differential equation (ODE) is a differential equation that contains one or more functions of one independent variable and its derivatives.

• The term ordinary is used in contrast with the term partial differential equation which may be with respect to more than one independent variable.

• The solutions to an ODE are functions that satisfy the equation.

• The order of an ODE is the number of the highest derivatives in the ODE

• Linear ODE: A linear ODE having independent variable t and the dependent variable y is an ODE of the form $a_0(t)y(n) + ... + a_{n-1}y' + a_n(t)y = f(t)$

• Examples of ODE:

• $(y'(x))^2 − x^2 = sin(y)$

• $y''(x) + 7y'(x)0 + 12y = e^x$

Solving ODE using SageMath

Example - Initial Value Problem (IVP)

• $x'(t) + x = 1$

• $x(0) = 2$

Sage's desolve may not be able to solve some difficult ODE.

Application: Falling Body

Consider a mass m falling due to gravity. We orient coordinates to that downward is positive. Let $x(t)$ denote the distance the mass has fallen at time t and $v(t)$ is the velocity at time t. Let assume only two forces act: the force due to gravity $F_{grav} = mg$, and the force due to air resistance $F_{res}$. So $F_{total} = F_{grav} + F_{res}$. We assume that air resistance is proportional to velocity: $F_{res} = -kv = -x'(t)$, where $k >= 0$ is a constant. By Newton's second law: $F_{total} = ma = mx''(t)$. Putting these all together gives $mx''(t) = mg - kx'(t)$, or $v'(t) + \frac{k}{m}v(t) = g$

Example: Suppose we have a mass 100kgs(with chute). The chute is released 30 secods after the jump from a height of 2000m. The force due to air resistence is giving by $F_{res} - -kv$, where $k=15$ when chute closed and $k=100$ when chute open. Suppose g = 9.8. Find

(a) The velocity functions during the time when the chute is closed (i.e $<= t <=30$ seconds)

The velocity at t = 30 is

(b) The velocity functions during the time when the chute is open (i.e., $30 <= t$ seconds)

Let's plot the results for both two periods

Numerical solutions - Euler's method

Euler's Method is intended for approximating solutions to differential equations that cannot be solved with a nice formula.

Say we have a generic initial value problem $\frac{dy}{dx} = f(x, y)$ with $y(x_0) = y_0$. We are not able to find a nice formula for the solution $y = g(x)$. We know if $f(x, y)$ is a continuous function, there exist a solution. Our goal is to say as much as we can about this solution, even though we may not be able to write down a formula for it.

Let's think about tangent line of the solution curve, since the tangent line is a good approximation to a curve near the point $(x_0, y_0)$. The equation of the tangent line to the solution at the point $(x_0, y_0)$ is $y=y_0 +f(x_0, y_0)(x-x_0)$. Let's move a short distance along the tangent line by increasing x by some small amount h, thwn we arrive at a new point, which we will call $(x_1, y_1)$. We can say the point $(x_1, y_1)$ is a almost on the solution curve $y = g(x)$. We can now find the tangent line to the solution at this new point $y = y_1 + f(x_1, y_1)(x-x_1)$, and then we increase x by the small amount h again, and come up with a new point $(x_2, y_2)$. Continuing with this method, we use the value of y found at each step to calculate the next value of y. This process can be stated by Euler's Formula: $y_{n+1} = y_n + f(x_n, y_n)(x_{n+1} - x_n)$

These codes use euler's method to find the numerical solution of the ODE $y'(x) = x^2-y+3$, with initial condition $y(0) = 1$, with step size 0.1 for x from 0 to 1

Plot of euler's method for the ODE $y'(x) = x^2 - y^3$, with intial condition $y(0) = 1$, with step size 1/3, 0.1, for x from 0 to 1

Remark: With smaller stepsize, the graph tends to be more smooth. The curve approximated will be more like the actual solution curve.

Homogeneous Systems of Linear ODE

Consider the system of differential equations

• $\frac{dx}{dt} = -y + 3x$

• $\frac{dy}{dt} = x - 5$

The general solution of this system

Now give some initial condition

scipy.integrate.ode

scipy.integrate.ode is a generic interface class to numeric integrators

Example:

• $\frac{dy}{dx} + y^2 = x^3, y(0)=1$

Exercise 1

Consider a tank with 200 liters of salt-water solution, 30 grams of which is salt. Pouring into the tank is a brine solution at a rate of 4 liters/minute and with a concentration of 1 grams per liter. The “well-mixed” solution pours out at a rate of 5 liters/minute. Find the amount at time t, and plot the amount for t in [0, 200]

You can learn more about the mixing problem in https://docs.google.com/document/d/1aLDQwp9zQ0KL90y7C8S5HbuYKDITEq2X_qUMK7w2Pkk/edit

Solution:

Exercise 2

Consider:

• $y'(x) = (y-1)(y+1)$

• $y(0) = 1/2$

Sage cannot solve this (implicit) solution for x(t). Show the implicit plot.

Solution:

Source

• https://www.usna.edu/Users/math/wdj/_files/documents/teach/sm212/DiffyQ/des-book-2009-11-24.pdf

• http://doc.sagemath.org/html/en/reference/calculus/sage/calculus/desolvers.html

• https://en.wikipedia.org/wiki/Ordinary_differential_equation