CODEWITHSUNDEEP
juliadifferential-equations
Sayeed
Sayeed

Reputation: 69

How to extract the evolution of the derivative of a variable from a differential equation solution in Julia

I am solving a differential equation in Julia the corresponding equations are given in the code below. Now by solving the differential equation what I am getting is the two variables only. But for my work, I also want the derivatives of the variables in each time step too after the integration is done, in which I am unable to proceed further.

using DelimitedFiles
using LightGraphs
using LinearAlgebra
using Random
using PyPlot
using BenchmarkTools
using SparseArrays

const N= 100;
Adj=readdlm("sf_simplicial_100.txt")
G=Graph(Adj)
global A = adjacency_matrix(G)
global deg=degree(G)
global omega=deg
mean(deg)


A2=zeros(N,N,N);
for i in 1:N
for j in 1:N
    for k in 1:N
        if (A[i,j]==1 && A[j,k]==1 && A[k,i]==1)
            A2[i,j,k]=1
            A2[i,k,j]=1
            A2[j,k,i]=1
            A2[j,i,k]=1
            A2[k,i,j]=1
            A2[k,j,i]=1
        end;
    end;
end;
end;

K= sum(p -> A2[:,:,p], 1:N)
deg_sim= sum(j -> K[:,j], 1:N)/2;
deg_sim2=2*deg_sim;

function kuramoto(du,u, pp, t)

u1 = @view u[1:N] #### θ
u2 = @view u[N+1:2*N]  ####### λ
du1 = @view du[1:N]  #### dθ
du2 = @view du[N+1:2*N]  ####### dλ


α1=0.08
β1=0.04
σ1=1.0
σ2=1.0
λ0=pp

####### local_order
z1 = Array{Complex{Float64},1}(undef, N)
mul!(z1, A, exp.((u1)im))
z1 = z1 ./ deg

####### generalized_local_order
z2 = Array{Complex{Float64},1}(undef, N)
z2= (diag(A*Diagonal(exp.((u1)im))*A*Diagonal(exp.((u1)im))*A))
z2 = z2 ./ deg_sim2

####### equ of motion
@. du1  = omega + u2 *( σ1 * deg * imag(z1 * exp((-1im) * u1)) + σ2 * deg_sim * imag(z2 * exp((-1im) * 2*u1)))
@. du2 = α1 *(λ0-u2)- β1 * (abs(z1)+ abs(z2))/2.0

return nothing
end;

using DifferentialEquations

# setting up time steps and integration intervals

 dt = 0.01 # time step
 dts = 0.1 # save time
 ti = 0.0
 tt = 1000.0 
 tf = 5000.0
 nt = Int(div(tt,dts))
 nf = Int(div(tf,dts))

 tspan = (ti, tf); # time interval
 pp=0.75
 ini=readdlm("N=100/initial_condition.txt")
 u0=[ini;pp*ones(N)];
 du = similar(u0);


 prob = ODEProblem(kuramoto,u0, tspan, pp)
 sol = solve(prob, RK4(), reltol=1e-4, saveat=dts,maxiters=1e10,progress=true)

Upvotes: 1

Views: 193

Answers (1)

Chris Rackauckas
Chris Rackauckas

Reputation: 19152

Use the derivative of the interpolation sol(t,Val{1}). To do this you'll want to not use saveat. Otherwise you can use the SavingCallback.

Upvotes: 1

Related Questions