# ODE boundary value problems

Let's have a look at how the DifferentialEquations.jl package handles boundary value problems. Let's solve the simple harmonic equation 
$$
u'' + 3 u = 0
$$
on the interval $t \in [0,2\pi]$ where $u(0)=7$ and $u(2 \pi)=0$. The exact solution is given by 
$$
\begin{align}
u(t) = 7 \cos \left(\sqrt{3}t \right) - 7 \cot \left(2\sqrt{3}\pi \right) \sin \left(\sqrt{3}t \right).
\end{align}
$$

The second order equation can be rewritten as the first order system
$$
\begin{align}
\left( u \right)' =& u', \\
\left( u' \right)' =& - 3 u.
\end{align}
$$
So let's solve this two point boundary value problem numerically using Julia.


```julia
#using Pkg
#Pkg.add( "BoundaryValueDiffEq" )
```


```julia
using DifferentialEquations, BoundaryValueDiffEq, Plots
using BenchmarkTools
```


```julia
function harmonic!(du,u,p,t)
  du[1] = u[2]
  du[2] = - 3 * u[1]
end
```




    harmonic! (generic function with 1 method)




```julia
function boundary_conditions!(residual, u, p, t)
    residual[1] = u[1][1] - 7
    residual[2] = u[end][1]
end
```




    boundary_conditions! (generic function with 1 method)




```julia
tspan = (0.0,2*pi)
tstep = 0.05
initial_guess = [7.0,-1.0]
bvp = TwoPointBVProblem( harmonic!, boundary_conditions!, initial_guess, tspan )
sol = solve( bvp, MIRK4(), dt=tstep);
```


```julia
plot( sol, vars=(1), lw=2, label="Numerical solution" )

function exact( t )
     return 7*( cos(sqrt(3)*t) - cot(2*sqrt(3)*pi) * sin(sqrt(3)*t) )
end

plot!( sol.t, t->exact(t), lw=2,ls=:dash, label="Exact solution" )
```




    

    




```julia
error = map( exact, sol.t )
for i in 1:length(sol.t)
    error[i] = error[i] - sol.u[i][1]
end
plot( sol.t, error, label="Error", lw=2 )
```




    

    



Let's try a slightly more difficult example such as the Blasius equation
$$
f''' + f f'' = 0,
$$
with the conditions $f(0)=f'(0)=0$ and $f'(\infty)=1$. We can rewrite the Blasius equation as the first order system
$$
\begin{align}
\left( f \right)' =& f', \\
\left( f' \right)' =& f'', \\
\left( f'' \right)' =& -f f''.
\end{align}
$$


```julia
function Blasius!(du,u,p,t)
    du[1] = u[2]
    du[2] = u[3]
    du[3] = -u[1] * u[3]
end
```




    Blasius! (generic function with 1 method)




```julia
function Blasius_bcs!(residual, u, p, t)
    residual[1] = u[1][1] 
    residual[2] = u[1][2]
    residual[3] = u[end][2] - 1
end
```




    Blasius_bcs! (generic function with 1 method)




```julia
tinf = 20
N = 200
tspan = (0.0,20.0)
tstep = (tspan[2] - tspan[1])/N
# Make an initial guess
t = tspan[1]:tstep:tspan[2]
function f_guess( t )
    return t * ( 1.0 - exp( -t ) )
end
function f_d_guess( t )
    return ( 1.0 - exp( -t ) ) + t * exp( -t )
end
function f_dd_guess( t )
    return 2 * exp( -t ) - t * exp( -t )
end
f = map( f_guess, t )
f_d = map( f_d_guess, t )
f_dd = map( f_dd_guess, t )
#initial_guess = [f,f_d,f_dd]
initial_guess = [1.0,1.0,0.0]
bvp_Blasius = TwoPointBVProblem( Blasius!, Blasius_bcs!, initial_guess, tspan )
@btime sol_Blasius = solve( bvp_Blasius, MIRK4(), dt=tstep );
```

      2.528 s (40327661 allocations: 3.56 GiB)



```julia
plot( sol_Blasius, xlims=(0,10), ylims=(0,2), label=["f" "f'" "f''"] )
```




    

    



Adding the initial guesses doesn't seem to help much with speed or the memory allocation. How can the performance be improved?


```julia

```

40aff6389cb506b54c13bcf290f9b97cfae8860f
ipynb
Jupyter Notebook
ODE_BVP.ipynb
anthonyoneill/jupyter-notebooks
9e1b0a6e1af51d38447adda01ac35060d1097c4e
MIT
ODE_BVP.ipynb
anthonyoneill/jupyter-notebooks
9e1b0a6e1af51d38447adda01ac35060d1097c4e
MIT
ODE_BVP.ipynb
anthonyoneill/jupyter-notebooks
9e1b0a6e1af51d38447adda01ac35060d1097c4e
MIT
Qwen/Qwen-72B
1. YES
2. YES