# Accuracy of TRBDF2 vs matlab's ode23tb?
**URL:** https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801
**Category:** Modelling & Simulations
**Tags:** differentialequation
**Created:** [October 2, 2024, 2:54am UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801 "2024-10-02T02:54:34Z")
**Posts on this page:** 16
**Page:** 1
### Author:  [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)
#### Post date: [October 2, 2024, 2:54am UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/1 "2024-10-02T02:54:34Z")
I noticed that trbdf2 solved my system the fastest but was much less accurate than matlab’s ode23tb. Some derivatives that should be zero are somehow non-zero even 100s of steps later.
By this I mean, `du[m:n]=0` from solution index=i to i+100 when calculated on the input u, but still `u[m:n]` keeps getting updated.
I am using TRBDF2 without auto-differentiation.
The issue goes away if I terminate the solver and restart the calculation when I am sure the derivatives wil be zero.
Is this a side effect of ‘smoothed derivatives’ mentioned in its documentation? When the derivatives are non-zero, they are pretty huge. Solvers like Rodas5P or Rosenbrock23 are correct but much slower. (Sometimes slower than using ode15 in matlab).
What other solvers should I try to use? My system isn’t big right now (~50 variables) but could be later.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [October 2, 2024, 6:00am UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/2 "2024-10-02T06:00:49Z")
Do you set the tolerance for each solvers?
---
### Author:  [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)
#### Post date: [October 2, 2024, 6:35am UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/3 "2024-10-02T06:35:16Z")
Yes each solver is supplied a reltol and abstol (of 1e-6)
---
### Author:  [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [October 2, 2024, 6:42am UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/4 "2024-10-02T06:42:11Z")
Can you provide some example code for the folks here to check?
---
### Author:  [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)
#### Post date: [October 2, 2024, 6:59am UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/5 "2024-10-02T06:59:27Z")
I’ll try and come up with a MWE, but it’ll take some time. I don’t know much about solvers, so I was wondering if someone encountered this before and had an opinion until I could make a mwe.
A naive pseudocode I have in mind could be something like
```julia
du[1:5] = rand(5)
du[6:8] = t in (10,12)? rand(3)*1e10 : 0
tspan = (0, 50)
```
But i’ll have to get back and run it
---
### Author:  [@danielwe](https://discourse.julialang.org/u/danielwe)
#### Post date: [October 2, 2024, 7:36am UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/6 "2024-10-02T07:36:21Z")
> [@pulk\_jain](#):
>
> ```julia
> du[1:5] = rand(5)
> du[6:8] = t in (10,12)? rand(3)*1e10 : 0
>
> ```
This is not an ODE. Is your actual system sampling random numbers?
---
### Author:  [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)
#### Post date: [October 2, 2024, 8:45am UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/7 "2024-10-02T08:45:33Z")
My actual system is not random numbers. I was just thinking of a system that’d recreate a regime with large derivatives and then zero derivatives later.
My system is more like two balloons floating in air and forces peak when they collide.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 2, 2024, 1:56pm UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/8 "2024-10-02T13:56:24Z")
This needs an example. With rands in there, I suspect there’s model code issues… hard to know what’s going on here without seeing code.
---
### Author:  [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)
#### Post date: [October 2, 2024, 8:33pm UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/9 "2024-10-02T20:33:54Z")
Okay, so I have a mwe example here. This model may not make much physical sense but I made this to recreate sharp gradients with small oscillation restricted to a certain domain. (My actual problem is much bigger and complicated than this)
Also, here the condition on u[3] is easily translatable to a condition in time as u[3] is just 3\*t.
```julia
function dynamics!(du,u,p,t)
du[1] = 2.0*u[3]
du[2] = 4.
du[3] = 3.
du[4] = ((u[3]>6) && (u[3]<20)) ? p*1e5*(u[2]-u[3])^2*u[1]*sin(t) : 0.0
du[5] = ((u[3]>6) && (u[3]<20)) ? 1e8*u[1]*cos(t) : 0.0
end
function condition(u,t,integrator)
u[3]>25.5 && norm(get_du(integrator)[4:5])!=0
end
function affect!(integrator)
@show integrator.t, integrator.u
terminate!(integrator)
end
cb = DiscreteCallback(condition,affect!)
prob = ODEProblem{true}(dynamics!, zeros(5), (0,10), 20);
sol = solve(prob, Tsit5(); callback = cb, dtmax=0.01, reltol =1e-6, abstol =1e-6, maxiters =2*1e5)
sol2 = solve(prob, TRBDF2(); callback = cb, dtmax=0.01, reltol =1e-6, abstol =1e-6, maxiters =2*1e5)
```
The Tsit5() solution is not stopped by the callback and runs smoothly while TRBDF2 is stopped. If TRBDF2 is run without callback, u[4] and u[5] still keep getting updated at t=10 (which is atleast 300 time steps away from u[3]=20).
With same tolerances matlab’s ode23tb does not run into these issues.
My question is if they are similar algorithms, there should not be such a big difference.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [October 3, 2024, 2:08am UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/10 "2024-10-03T02:08:45Z")
This seems to be plenty fast and working just fine as long as you don’t expect the du to be exactly 0, but rather a number \ @btime sol = solve(prob, Rodas5P(); reltol =1e-6, abstol =1e-6);
432.221 μs (674 allocations: 66.91 KiB)
julia> @btime sol = solve(prob, TRBDF2(); reltol =1e-6, abstol =1e-6);
381.361 μs (580 allocations: 53.39 KiB)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 5, 2024, 3:33pm UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/11 "2024-10-05T15:33:17Z")
Is this solved? The issue was just that `dtmax` was set really low. I presume you didn’t do that on the MATLAB side? I presume that instead of `dtmax` you instead intended to use `saveat`.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [October 5, 2024, 4:09pm UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/12 "2024-10-05T16:09:07Z")
I believe it’s solved. The one thing that is slightly interesting is that `get_du` for `TRBDF2` is a bit inaccurate for the last 2 equations. I think that might be an inherent property of the the newton solver (i.e. it stops improving once it’s found a solution within tolerance) but it is somewhat surprising that when your equation is `du[5] = 0` `getdu` doesn’t return `0` exactly.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 5, 2024, 4:18pm UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/13 "2024-10-05T16:18:08Z")
I think it’s partially because this method does not have a great interpolation on it, and we should improve the interpolation choice here.
---
### Author:  [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)
#### Post date: [October 7, 2024, 8:44pm UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/14 "2024-10-07T20:44:16Z")
I’d say its solved if non-zero gradients for last two equations is expected behavior for TRBDF2. Thanks @Oscar_Smith for your inputs.
My main issue stemmed from expecting exact zero derivatives (because I intend to use that as a callback condition later), for which I’ve switched to non-stiff solvers just to be sure of the accuracy. I’ll then again try TRBDF2 and compare derivatives to a tolerance rather than 0, to see if the results match.
> I presume you didn’t do that on the MATLAB side? I presume that instead of `dtmax` you instead intended to use `saveat` .
I use the same dtmax on both platforms. I want to set a dtmax to resolve the timescales associated with the problem. I believe that must have been the reason to have that feature in the first place?
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [October 7, 2024, 10:26pm UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/15 "2024-10-07T22:26:52Z")
`dtmax` should almost never be used (it’s only use is when debugging/working with solvers with buggy error control). Unless you have very strong reasons to distrust the accuracy of your solver, you should not set `dtmax` and should let the adaptivity figure out appropriate dt on it’s own.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 8, 2024, 10:17am UTC](https://discourse.julialang.org/t/accuracy-of-trbdf2-vs-matlabs-ode23tb/120801/16 "2024-10-08T10:17:52Z")
> [@pulk\_jain](#):
>
> My main issue stemmed from expecting exact zero derivatives (because I intend to use that as a callback condition later), for which I’ve switched to non-stiff solvers just to be sure of the accuracy. I’ll then again try TRBDF2 and compare derivatives to a tolerance rather than 0, to see if the results match.
Rosenbrock methods should be able to enforce exact zeros better than an implicit method.