# Imposing constraints on solutions by DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/imposing-constraints-on-solutions-by-differentialequations-jl/104255>\
**Category:** Numerics\
**Tags:** differentialequation\
**Created:** [September 26, 2023, 11:03am UTC](https://discourse.julialang.org/t/imposing-constraints-on-solutions-by-differentialequations-jl/104255 "2023-09-26T11:03:50Z")\
**Posts on this page:** 10\
**Page:** 2

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 27, 2023, 12:34pm UTC](https://discourse.julialang.org/t/imposing-constraints-on-solutions-by-differentialequations-jl/104255/22 "2023-09-27T12:34:18Z")

</div>

I don’t get what you’re saying. What do you mean by “doesn’t work” and “these numbers”?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [September 27, 2023, 12:38pm UTC](https://discourse.julialang.org/t/imposing-constraints-on-solutions-by-differentialequations-jl/104255/23 "2023-09-27T12:38:22Z")

</div>

if it is linear, why not try a exponential method?

---

<div class="post-metadata">

**Author:** ![NitinGupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nitingupta/32/53176_2.png) [@NitinGupta](https://discourse.julialang.org/u/NitinGupta)\
**Post date:** [September 27, 2023, 1:01pm UTC](https://discourse.julialang.org/t/imposing-constraints-on-solutions-by-differentialequations-jl/104255/24 "2023-09-27T13:01:21Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![NitinGupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nitingupta/32/53176_2.png) [@NitinGupta](https://discourse.julialang.org/u/NitinGupta)\
**Post date:** [September 27, 2023, 1:02pm UTC](https://discourse.julialang.org/t/imposing-constraints-on-solutions-by-differentialequations-jl/104255/25 "2023-09-27T13:02:26Z")

</div>

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`.

---

<div class="post-metadata">

**Author:** ![NitinGupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nitingupta/32/53176_2.png) [@NitinGupta](https://discourse.julialang.org/u/NitinGupta)\
**Post date:** [September 27, 2023, 1:15pm UTC](https://discourse.julialang.org/t/imposing-constraints-on-solutions-by-differentialequations-jl/104255/26 "2023-09-27T13:15:18Z")

</div>

here’s a MWE -

```julia
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!

 ![Screenshot 2023-09-27 at 3.10.49 PM](https://global.discourse-cdn.com/julialang/original/3X/c/4/c46ccbf88a24eec00d382532be460dd22ab96f37.png)  
 ![Screenshot 2023-09-27 at 3.14.26 PM](https://global.discourse-cdn.com/julialang/original/3X/7/4/741fd432688085e55dfa3e84f524504c3338c235.png)

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [September 27, 2023, 1:26pm UTC](https://discourse.julialang.org/t/imposing-constraints-on-solutions-by-differentialequations-jl/104255/27 "2023-09-27T13:26:35Z")

</div>

> [@NitinGupta](#):
>
> `dΦ[MaxCoeff] = -im*(real(an[MaxCoeff-1]) * Φ[MaxCoeff] + real(bn[MaxCoeff-1]) * Φ[MaxCoeff-1] + real(bn[MaxCoeff]) * Φ[MaxCoeff+1])`

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

---

<div class="post-metadata">

**Author:** ![NitinGupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nitingupta/32/53176_2.png) [@NitinGupta](https://discourse.julialang.org/u/NitinGupta)\
**Post date:** [September 27, 2023, 1:33pm UTC](https://discourse.julialang.org/t/imposing-constraints-on-solutions-by-differentialequations-jl/104255/28 "2023-09-27T13:33:04Z")

</div>

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 term`bn[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.

---

<div class="post-metadata">

**Author:** ![NitinGupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nitingupta/32/53176_2.png) [@NitinGupta](https://discourse.julialang.org/u/NitinGupta)\
**Post date:** [September 27, 2023, 1:34pm UTC](https://discourse.julialang.org/t/imposing-constraints-on-solutions-by-differentialequations-jl/104255/29 "2023-09-27T13:34:19Z")

</div>

for example, if I use the following values (setting `MaxCoeff=10`) then I get huge errors :

```julia
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]

```

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [September 27, 2023, 1:40pm UTC](https://discourse.julialang.org/t/imposing-constraints-on-solutions-by-differentialequations-jl/104255/30 "2023-09-27T13:40:04Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![NitinGupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nitingupta/32/53176_2.png) [@NitinGupta](https://discourse.julialang.org/u/NitinGupta)\
**Post date:** [September 27, 2023, 2:56pm UTC](https://discourse.julialang.org/t/imposing-constraints-on-solutions-by-differentialequations-jl/104255/31 "2023-09-27T14:56:39Z")

</div>

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 🙂

[Previous page](https://discourse.julialang.org/t/imposing-constraints-on-solutions-by-differentialequations-jl/104255.md?page=1)