Free_ion
Free_ion

Reputation: 97

Defining a function with the solution of a differential equation in Julia

I want to use the solution of a differential equation into a piecewise function, but when plotting the result Julia returns the error MethodError: no method matching Float64(::LinearAlgebra.Transpose{Float64, Vector{Float64}}). With the code below, you can see the plot works fine until r=50, and that the error shows up for r>50. How can I use the result of the differential equation in a function?

using Plots, DifferentialEquations

const global Ωₘ = 0.23
const global Ωᵣ = 0.0001
const global Ωl = 0.67
const global H₀ = 70.

function lum(du,u,p,z)
    du[1] = -u[1]/(1+z) + ((1/H(z))-1/H₀)/(1+z)
end
da0 = [0.]
zspan = (0.,10.)
lumprob = ODEProblem(lum,da0,zspan)
lumsol = solve(lumprob)

α1(r) = π- asin(3*sqrt(3)*sqrt(1-2/r)/r)
α2(r) = asin(3*sqrt(3)*sqrt(1-2/r)/r)
α3(r) = 3*sqrt(3)/r
α4(r) = 3*sqrt(3)/lumsol(r)

function αtot(r)
    if 1.99<r<3.0
        α1(r)
    elseif 3.0<r<=10.0
        α2(r)
    elseif 10.0<r<=50.0
        α3(r)
    elseif 50.0<r
         α4(r)
    end
end
p2=plot(αtot,xlims=(2.0,51.0),xaxis=:log,xlabel="Rₒ/m",label="α")

Upvotes: 2

Views: 181

Answers (1)

Vitaliy Yakovchuk
Vitaliy Yakovchuk

Reputation: 493

In your case lumsol returns a 1-element vector but not a number, thus you need to change line

α4(r) = 3*sqrt(3)/lumsol(r)

to the

α4(r) = 3*sqrt(3)/lumsol(r)[1]

(or refactor the code to take into account that fact)

The function plot should take a function that takes a number and returns a number for the case, your αtot function returns a vector if 50 < r.

Upvotes: 3

Related Questions