# ODE solving

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.

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

`T1` isn’t defined in your example.

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

Hi,

I added the temperatures from the text.

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

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

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.

Here is the scipy variant:

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

Well that’s easy to check:

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

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:

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

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

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

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

1 Like