# MAE 103: Differential Equations in Julia

## Solve differential equations in Julia and pathlines

In today's lecture, we will learn how to solve ordinary differential equations. We will then use it to plot pathlines of fluid particles

### A 1D differential equation

In this first example, we will solve the equation: $\frac{d u}{dt} = -u$ with $u(t=0) = 1$

We know the analytical solution for this problem : $u(t) = exp(-t)$

In Julia, we can call the package OrdinaryDiffEq to solve numerically ODEs

using DifferentialEquations
Shift+Enter to run

The first step is to define the right hand side of the differential equation, here f(u,p,t) = -u

f(u,p,t) = -u
Shift+Enter to run
nil

The initial condition is u0 = 1

u0=1.0
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
Shift+Enter to run
nil

The variable prob holds all the useful informations about the differential problem : right-hand side, initial condition, time span...

We will solve this ODe with a Runge-Kutta 4 scheme

sol = solve(prob,RK4())
Shift+Enter to run
nil

#### Plotting the solution

using Plots
Shift+Enter to run
plot(sol.t, exp.(-sol.t), linewidth = 4, label = "Analytical solution")
plot!(sol,linewidth=4, xaxis="t",yaxis="u", linestyle =:dash, label = "Numerical solution" ) # legend=false
Shift+Enter to run
image/svg+xml
using Interpolations
Shift+Enter to run

### Solve for the pathline of a fluid particle in a Couette Flow

For a Couette flow, the velocity profile is: $\boldsymbol{u} = u(y)\boldsymbol{e}_x = U\frac{y}{H} \boldsymbol{e}_x$ with $U$ the velocity of the upper plate and $H$ the depth of the gap

The trajectory of a fluid particle A is denoted $\boldsymbol{x}_A$ and a solution of:

with initial condition $\boldsymbol{x}_A(t= 0) = \boldsymbol{x}_0$

Let's solve for the pathline of a fluid particle. For this problem we have :

Define the parameters for the Couette flow

U = 1
H = 1
couette_velocity(y) = U*y/H
Shift+Enter to run
nil

Let's define the right-hand side of the equation

function couette(du,u,p,t)
 du[1] = couette_velocity(u[2]) # U * (yₐ/H)
 du[2] = 0.0
end
Shift+Enter to run
nil
## Initial position of the fluid particle
x0=[0.0, 0.5]
# Time span
tspan = (0.0,1.0)
prob_couette = ODEProblem(couette, x0, tspan)
Shift+Enter to run
nil

Now, we can solve the differential equation and get the pathline of this fluid particle

couette_sol  = solve(prob_couette, RK4())
Shift+Enter to run
nil
plot(couette_sol, vars = 1, legend = false, xlabel = "t", ylabel = "xA(t)")
Shift+Enter to run
image/svg+xml
plot(couette_sol, vars = 2, legend = false, xlabel = "t", ylabel = "yA(t)")
Shift+Enter to run
image/svg+xml
plot(couette_sol, vars = (1,2), ylim =(0, H))
plot!([1.0], line = :hline, linewidth = 8, label = "Upper wall")
plot!([0.0], line = :hline, linewidth = 8, label = "Lower wall")
scatter!([x0[1]], [x0[2]])
Shift+Enter to run
image/svg+xml

## Trajectory of a fluid particle in a cellular flow

A = 32.0
B = 3.2
ω = 9.2
ψ(x, y, t) = A*cos(x + B*sin(ω*t))*cos(y)
ucell(x, y, t) =  A*cos(x + B*sin(ω*t))*sin(y)
vcell(x, y, t) = -A*sin(x + B*sin(ω*t))*cos(y)
Shift+Enter to run
nil
xrange = -4.0:0.05:12.0
yrange = -4.0:0.05:12.0
function ψtab(t)
    out  = zeros(xrange.len ,yrange.len)
for (i, xi) in enumerate(xrange)
    for (j, yj) in enumerate(yrange)

        out[i, j] = ψ(xi, yj , t)

    end
end

    return out
end
Shift+Enter to run
nil
plot(xrange, yrange, ψtab(0.0)', ratio = 1, xlim = (xrange[1], xrange[end]), ylim = (yrange[1], yrange[end]), colorbar = false,
            xlabel ="x", ylabel = "y", title = "Streamlines at t = "*string(0.0))
Shift+Enter to run
image/svg+xml
plot(xrange, yrange, ψtab(1.0)', ratio = 1, xlim = (xrange[1], xrange[end]), ylim = (yrange[1], yrange[end]), colorbar = false,
            xlabel ="x", ylabel = "y", title = "Streamlines at t = "*string(1.0))
Shift+Enter to run
image/svg+xml
plot(xrange, yrange, ψtab(2.0)', ratio = 1, xlim = (xrange[1], xrange[end]), ylim = (yrange[1], yrange[end]), colorbar = false,
            xlabel ="x", ylabel = "y", title = "Streamlines at t = "*string(2.0))
Shift+Enter to run
image/svg+xml
function cellular_flow(du,u,p,t)
 du[1] = ucell(u[1], u[2], t)
 du[2] = vcell(u[1], u[2], t)
end
Shift+Enter to run
nil
## Initial position of the fluid particle
x0_cell=[3.20, 1.88]
# Time span
tspan = (0.0,2.0)
prob_cellular = ODEProblem(cellular_flow, x0_cell, tspan)
Shift+Enter to run
nil
cellular_sol  = solve(prob_cellular, RK4())
Shift+Enter to run
nil
plot(xrange, yrange, ψtab(0.0)', ratio = 1, xlim = (xrange[1], xrange[end]), ylim = (0.0, 5.0), colorbar = false,
            xlabel ="x", ylabel = "y", title = "Streamlines at t = "*string(0.0))
plot!(cellular_sol, vars = (1,2))
scatter!([x0_cell[1]], [x0_cell[2]])
Shift+Enter to run
image/svg+xml
Shift+Enter to run
Shift+Enter to run