# ODE solving

**URL:** <https://discourse.julialang.org/t/ode-solving/99520>\
**Category:** Modelling & Simulations\
**Tags:** question, ode, nlsolve\
**Created:** [May 28, 2023, 4:52pm UTC](https://discourse.julialang.org/t/ode-solving/99520 "2023-05-28T16:52:32Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [May 28, 2023, 4:52pm UTC](https://discourse.julialang.org/t/ode-solving/99520/1 "2023-05-28T16:52:32Z")

</div>

I am Trying to model some temperature curves with a physically informed NN (PINN).  
I started solve the FOPDT first, to be able to apply it later.

```julia

function FOPDT(z)
    τ, Κ, θ = z
    return [T1[i] .- (1 ./ τ) .* (-T1[i] .+ Κ .* 100 .*(i .- θ)) for i in [1, 75, 130]]

end
initial_x = [1.0, 1.0, 1.0]
sol = nlsolve(FOPDT, initial_x, xtol=1e-2, ftol=1e-2, iterations=1_000_000)
sol

```

The temperatures at the time steps are:

0 s → 21.5 °C  
75 s → 46.0 °C  
130 s → 60.3 °C

But I was not able to find a static solution. Hence I used Scipy.optimize.fsolve. I tried to use  
this solution as starting parameters but this also did not converge. Are there some other hints, how to solve the equation? Did I do any other mistake, e.g. in function?

---

<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:** [May 28, 2023, 9:20pm UTC](https://discourse.julialang.org/t/ode-solving/99520/2 "2023-05-28T21:20:18Z")

</div>

`T1` isn’t defined in your example.

```julia
using NonlinearSolve

function FOPDT(dz,z,p)
  τ, Κ, θ = z
  return dz .= [T1[i] .- (1 ./ τ) .* (-T1[i] .+ Κ .* 100 .*(i .- θ)) for i in [1, 75, 130]]
end
initial_x = [1.0, 1.0, 1.0]
prob = NonlinearProblem(FOPDT, initial_x)
sol = solve(prob, NonlinearSolve.NewtonRaphson(), abstol=1e-2)

```

Can you please give enough information to run this?

---

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [May 28, 2023, 10:26pm UTC](https://discourse.julialang.org/t/ode-solving/99520/3 "2023-05-28T22:26:30Z")

</div>

Hi,

I added the temperatures from the text.

```julia

function FOPDT(z)
    τ, Κ, θ = z
   return [T - (1 / τ) * (-T + Κ * 100 *(i - θ)) for (T, i) in zip([21.543, 46.035,60.279], [1, 75, 130])]

end

initial_x = [145.7, 66.3, 73.4] # from scipy.optimize.fsolve
sol = nlsolve(FOPDT, initial_x, xtol=1e-1, ftol=1e-2)

```

The result is:

```julia
Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [145.7, 66.3, 73.4]
 * Zero: [-7.290966107926511e11, -2.2018393026492186e9, -72.45793039799061]
 * Inf-norm of residuals: 1.503362
 * Iterations: 183
 * Convergence: true
   * |x - x'| < 1.0e-01: true
   * |f(x)| < 1.0e-02: false
 * Function Calls (f): 184
 * Jacobian Calls (df/dx): 100

```

---

<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:** [May 28, 2023, 10:31pm UTC](https://discourse.julialang.org/t/ode-solving/99520/4 "2023-05-28T22:31:39Z")

</div>

That point very much evaluates to non-zero. But can you share the SciPy code? It would take the comparison in order to see what’s different in the function.

---

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [May 28, 2023, 11:00pm UTC](https://discourse.julialang.org/t/ode-solving/99520/5 "2023-05-28T23:00:15Z")

</div>

Here is the scipy variant:

```python
from scipy.optimize import fsolve

def f(z):
    tau, kappa, theta = z
    return [T-1/tau*(-T + kappa*100*(i - theta)) for T, i in zip([21.543, 46.035, 60.279], [0, 75, 130])]
tau, kappa, theta = fsolve(f, [1, 1, 1], xtol=1e-01, maxfev=100_000)

```

A solution was only found with adding the xtol parameter and setting it to: 1e-1. I don’t know how much to trust the result.

---

<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:** [May 29, 2023, 2:31am UTC](https://discourse.julialang.org/t/ode-solving/99520/6 "2023-05-29T02:31:43Z")

</div>

> [@Ordens\_Ritter](#):
>
> I don’t know how much to trust the result.

Well that’s easy to check:

```julia
>>> def f(z):
... tau, kappa, theta = z
... return [T-1/tau*(-T + kappa*100*(i - theta)) for T, i in zip([21.543, 46.035, 60.279], [0, 75, 130])]
... 
>>> tau, kappa, theta = fsolve(f, [1, 1, 1], xtol=1e-01, maxfev=100_000)
>>> f([tau, kappa, theta])
[3360.9606761659174, -27.294401084985765, -2515.757097939305]

```

It’s also not a zero in the Python code.

---

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [May 29, 2023, 2:12pm UTC](https://discourse.julialang.org/t/ode-solving/99520/7 "2023-05-29T14:12:53Z")

</div>

Strange. But I assume I did a mistake in the formula.  
In addition, I solved the solution for cooling phase and than for the heating phase.

Cooling:

```python
def f_off(z):
    _tau = z[0]
    def k(c: float):
        return c + 273.15

    return [dT + k(T)/_tau for dT, T in zip([-0.2425,], [76.006,] )]
r = fsolve(f_off, [1.0,], xtol=1e-6, maxfev=1000_000)

```

Heating:

```python
def f_on(z):
    _tau = 1439.818557
    _kappa, _theta = z
    def k(c: float):
        return c + 273.15

    return [dT + 1/_tau*(T + _kappa*100*(_theta -i)) for dT, T, i in zip([0.279, 0.3065], [46.035, 60.279], [75, 130])]
kappa, theta = fsolve(f_on, [66, 73], xtol=5e-00, maxfev=1000_000)

```

I got:  
tau = 1439.818557,  
kapp = 0.0097885  
theta = 72.93

Again no solution found in the Julia. Should I try another library? Which?

---

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [May 30, 2023, 8:48pm UTC](https://discourse.julialang.org/t/ode-solving/99520/8 "2023-05-30T20:48:03Z")

</div>

I did a bug fix and compared the python and the julia solution.  
Both are stable and results are the same with:

```julia
function FOPDT_off(z)
    τ = z[1]
	k(c) = c + 273.15
	
    return [dT + k(T) / τ for (dT, T) in zip([-0.2425,], [76.006,])]
end

```

```julia
function FOPDT_on(z)
        τ = 1439.8
	Κ, θ = z
	k(c) = c + 273.15
	
    return [dT + (1 / τ) * (k(T) + Κ * 100 *(θ -i)) for (dT, T, i) in zip(
		[0.279, 0.3065],
		[46.035,60.279], 
		[75, 130])]

end

```

The solutions are: τ=1439.8, Κ=0.009788, θ=-661.4 .