Help to adapt the resolution of an ODE in Matlab to Julia

Specific Domains Numerics
1

Hello guys. I’m having trouble adapting a Matlab code that resolves a second-order ODE to the Julia equivalent.

This is the code in Matlab. It plots a graph with 3 curves:

syms y(x);
dy = diff(y);
edo = diff(y,x,2) == -0.1 * dy - 0.8 * y;
Bound_1 = y(0) == 5e-3;
Bound_2 = dy(0) == 0;
conds = [Bound_1 Bound_2];
t = [0, 30];
ySol(x) = dsolve(edo, conds);
figure(1)
fplot(ySol, t)
hold on
fplot(diff(ySol), t)
fplot(diff(ySol,2), t)
hold off

Captura de tela de 2022-09-16 19-17-25
Captura de tela de 2022-09-16 19-17-251659×1376 274 KB

This is my equivalent Julia version. However, the resulting graph has only two curves.

f2(du, u, p, t) = -0.1du - 0.8u
u₀ = 5e-3
du₀ = 0.0
t = (0.0, 30.0)
prob = SecondOrderODEProblem(f2, du₀, u₀, t)
sol = solve(prob)
plot(sol)

image
image1244×810 44.9 KB

How do I make the u'' curve appear in the graph generated by the julia?

2

If you decrease the tolerances do they become the same?

3

The graph is correct but the curve representing the second order derivative is missing.

2 Likes
4

To interpolate the second-order derivative, you can use this:

ts = LinRange(t[1], t[2], 200)
z = sol(ts, Val{0})
dz = sol(ts, Val{1})

using Plots 
plot(sol)
plot!(ts, dz[1,:])

# or 
plot(ts, [z[2,:] z[1,:] dz[1,:]] )

The sol(...) syntax is used for interpolation, see Solution Handling · DifferentialEquations.jl
and the Val{1} indicates that we want the first derivative of the solution.

Since you have a second-order problem, the solution vector has the form z = (u', u) and accordingly we have z' = (u'', u'). Therefore we get dz[1,:] ~ u''.

2 Likes