I don’t get what you’re saying. What do you mean by “doesn’t work” and “these numbers”?
if it is linear, why not try a exponential method?
1 Like
By “doesn’t work” I mean that I get a large amount of error and by “these numbers” I refer to an and bn (from the image in an earlier comment) - basically, if I remove these numbers from the problem then I get really low error. So I think the issue must be with these numbers and I see nothing wrong with them!
you mean something like LawsonEuler? I tried several of them but they give Cannot create a dual over scalar type ComplexF64 even with autodiff=false.
here’s a MWE -
using Plots, OffsetArrays, LaTeXStrings
using DifferentialEquations
MaxCoeff=30
an=ones(MaxCoeff+1)
bn=ones(MaxCoeff)
an=OffsetVector(an,0:MaxCoeff)
function ComplexityWFSolver!(dΦ, Φ, p, t)
dΦ[1] = -im*(real(an[0]) * Φ[1] + real(bn[1]) * Φ[2])
dΦ[MaxCoeff] = -im*(real(an[MaxCoeff-1]) * Φ[MaxCoeff] + real(bn[MaxCoeff-1]) * Φ[MaxCoeff-1] + real(bn[MaxCoeff]) * Φ[MaxCoeff+1])
dΦ[MaxCoeff+1] = -im*(real(an[MaxCoeff]) * Φ[MaxCoeff+1] + real(bn[MaxCoeff]) * Φ[MaxCoeff])
for i in 2:MaxCoeff-1
dΦ[i] = -im*(real(an[i-1]) * Φ[i] + real(bn[i+1]) * Φ[i+1] + real(bn[i]) * Φ[i-1])
end
end
Φ₀ = zeros(Complex{Float64},MaxCoeff+1)
Φ₀[1] += 1.0
Δt = (0.0,10.0);
prob = ODEProblem(ComplexityWFSolver!, Φ₀, Δt, nothing);
sol = solve(prob,Tsit5(),reltol=1e-6,abstol=1e-6)
timelist=sol.t
sizetimelist=size(timelist)[1]
ComplexityWFProb = zeros(Float64,sizetimelist)
for i in 1:sizetimelist
for j in 1:MaxCoeff
ComplexityWFProb[i] += (abs(sol(timelist[i])[j])^2)
end
end
ProbabilityPlot=plot((timelist,ComplexityWFProb),
marker = (:circle), markersize= 0.75, label=L"${\sum_n ~ |ϕ_n|^2}$",
xlabel = L"t", ylabel = L"${\sum_n ~ |ϕ_n|^2}$")
display(ProbabilityPlot)
This produces the following, which is excellent because the y-axis should remain around 1 at all times. But now, instead of an=ones(MaxCoeff+1) and bn=ones(MaxCoeff) I use other values (they are BigFloat, I obtain them from solving my physical problem and I’m 100% certain they are correct because other people have also independently got the same values using different methods and different programming languages), I get the second picture which is terrible!
should be in the loop over the interior points and seems actually inconsistent with the main loop on coefficients. I think there is a shift in indexes, which would explain that it works perfectly when all coefficients are the same
Actually, it’s supposed to be outside because of the bounds on the vector bn because if you run the loop until i=MaxCoeff then the second termbn[i+1] will throw error. It’s consistent due to the way problem is constructed (believe me on this). It’s easy to see from the picture I posted in one of the earlier comments where the system is represented as a matrix equation.
for example, if I use the following values (setting MaxCoeff=10) then I get huge errors :
a = [0.75,2.25,
2.7930689840936798,
4.978292362076318,
4.852205732976179,
7.312243143858942,
7.734831978392002,
8.718418701931046,
11.246470922104564,
10.144237373425446,
13.839111980530182]
b = [0.7516648189186456,
2.498459702487865,
3.4790750299364706,
4.928722354966566,
6.433152863922314,
7.164134095776164,
9.383715606024372,
9.719953246001241,
11.697936414921527,
12.948876534350841]
There are only two special lines in your system, the first and the last which demand a special treatment.
All other should fit into the same loop and if there is a problem with an index in that case, it probably means you got an index wrong
1 Like
I’m stupid -_-. I made an index error and correcting it fixes the issue I was having. Although without implementing what the title of my question says. So I’ll accept Chris’s comment about using symplectic integrators because of the title.
Huge thanks to JM_Beckers and all others for patiently helping me 
1 Like