I have a system with a dependent mass matrix. It always return me some mensage of instability. Is it an intrisic characteristic of the system?
Any thing I can do about?
using DifferentialEquations
alpha=0.105
N0 = 5
Ar=2
Eeq=(1.0+N0)/(2.0*N0)+((9.0+6.0*N0+9.0*N0^2.0)^0.5/6.0/N0)
function update_func(A,u,p,t)
E = (9*N0-(u[1]+u[2]+u[3]))/(9*N0-3*(u[1]+u[2]+u[3]))
De = 0.1
A[1,1] = De*u[5]
A[1,2] = 0
A[1,3] = 0
A[1,4] = 0
A[1,5] = -2*De*u[1]
A[1,6] = 0
A[2,1] = 0
A[2,2] = De*u[5]
A[2,3] = 0
A[2,4] = -2*De*u[1]*u[5]/u[4]
A[2,5] = 0
A[2,6] = 0
A[3,1] = 0
A[3,2] = 0
A[3,3] = De*u[5]
A[3,4] = 0
A[3,5] = 0
A[3,6] = -2*De*u[3]*u[5]/u[6]
A[4,1] = 0
A[4,2] = 0
A[4,3] = 0
A[4,4] = 1.
A[4,5] = 0
A[4,6] = 0
A[5,1] = 0
A[5,2] = 0
A[5,3] = 0
A[5,4] = 1/u[4]
A[5,5] = 1/u[5]
A[5,6] = 1/u[6]
A[6,1] = 1/(u[1]-u[3]) + 2*N0/E/(3*N0-u[1]-u[2]-u[3])^2
A[6,2] = 2*N0/E/(3*N0-u[1]-u[2]-u[3])^2
A[6,3] = -1/(u[1]-u[3]) + 2*N0/E/(3*N0-u[1]-u[2]-u[3])^2
A[6,4] = 1/u[4]
A[6,5] = 0
A[6,6] = 1/u[6]
end
function casting(du,u,p,t)
E = (9*N0-(u[1]+u[2]+u[3]))/(9*N0-3*(u[1]+u[2]+u[3]))
kx = ((1.0-alpha)+alpha*E*u[1])*(1-E*u[1])
ky = ((1.0-alpha)+alpha*E*u[2])*(1-E*u[2])
kz = ((1.0-alpha)+alpha*E*u[3])*(1-E*u[3])
du[1] = kx
du[2] = ky
du[3] = kz
du[4] = -Ar*((u[2]-u[3])/(u[1]-u[3]))^(1/2)
du[5] = 0
du[6] = 0
nothing
end
guess = 2
uu0 = [guess; 1/Eeq; 1/Eeq; 1.0; 1.0; 1.0]
M = DiffEqArrayOperator(ones(6,6),update_func=update_func)
tspan = (0, 1.0)
prob = ODEProblem(ODEFunction(casting, mass_matrix=M), uu0, tspan)
sol = solve(prob,Rodas5())
┌ Warning: Instability detected. Aborting
└ @ SciMLBase C:\Users\jorgem63\.julia\packages\SciMLBase\n3U0M\src\integrator_interface.jl:351
retcode: Unstable
Interpolation: 3rd order Hermite
t: 1-element Vector{Float64}:
0.0
u: 1-element Vector{Vector{Float64}}:
[2.0, 0.8759615953640396, 0.8759615953640396, 1.0, 1.0, 1.0]