# Dependent mass matrix with RadauIIA5

**URL:** <https://discourse.julialang.org/t/dependent-mass-matrix-with-radauiia5/43189>\
**Category:** General Usage\
**Tags:** diffeq\
**Created:** [July 16, 2020, 4:34pm UTC](https://discourse.julialang.org/t/dependent-mass-matrix-with-radauiia5/43189 "2020-07-16T16:34:20Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jorge\_Silva](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jorge_silva/32/14976_2.png) [@Jorge\_Silva](https://discourse.julialang.org/u/Jorge_Silva)\
**Post date:** [July 16, 2020, 4:34pm UTC](https://discourse.julialang.org/t/dependent-mass-matrix-with-radauiia5/43189/1 "2020-07-16T16:34:20Z")

</div>

When running this code Im obtaing not expect not smooth solutions.

> ```
> using Plots
> using DifferentialEquations
> 
> Ta=20.0+273.15;                                                               
> ka=1.5207e-11*Ta^3 - 4.8574e-8*Ta^2 + 1.0184e-4*Ta - 3.9333e-4;                                     
> ro_a=360.77819/Ta; 
> eta_a=(-1.1555e-14*Ta^3 + 9.5728e-11*Ta^2 + 3.7604e-8*Ta - > 3.4484e-6)*ro_a;   
> Q=0.6/3600;                
> R0=0.00135;                
> L=0.140;                   
> T0=150.0+273.15;            
> ro=781.5;                        
> N0=3.0;                     
> eta_T0=2.2e4;              
> lambda_T0=2.0;             
> Ea=5.77e4;                
> alpha=0.105;               
> g=9.81;                     
> u0=Q/(ro*3.1415*R0^2);   
> G0=eta_T0/lambda_T0; 
> V1=0.42*L*ka/(ro*u0^(2/3)*R0^(5/3))*(ro_a/eta_a)^(1/3); 
> V3=G0/(ro*T0); 
> C1=ro*u0^2/G0;
> C3=(ro*g*L)/G0; 
> Eeq=(1.0+N0)/(2.0*N0)+((9.0+6.0*N0+9.0*N0^2.0)^0.5/6.0/N0);
> czz0=2.0;
> 
> function filament(du,u,p,t)
> 
> Cp=1980+3.7*T0*u[5] 
> EE=(9.0*N0-2.0*u[1]-u[2])/(9.0*N0-6.0*u[1]-3.0*u[2])
> dEE=18.0*N0/(9.0*N0-6.0*u[1]-3.0*u[2])^2
> lambda=lambda_T0*exp(Ea/(8.314*T0)*((1.0/u[5])-1.0))
> Cf=0.205*((ro_a*u[4]*u0*u[3]*R0)/eta_a)^(-0.61)
> C2=2.0*Cf*ro_a*L*u0^2/(G0*R0)   
> 
> println(u[3])
> 
> du[1] = u[5]/(lambda*u0/L)*((1.0-alpha) + alpha*EE*u[1]/u[5])*(EE*u[1]/u[5]-1.0)
> du[2] = u[5]/(lambda*u0/L)*((1.0-alpha) + alpha*EE*u[2]/u[5])*(EE*u[2]/u[5]-1.0)
> du[3] = 0.0
> du[4] = -C2*u[4]^2/u[3] + C3
> du[5] = -V1/(Cp*u[4]^(2/3)*u[3]^(5/3))*(u[5]-Ta/T0)
> 
> nothing
> end    
> 
> function update_func(M,u,p,t)
> 
> Cp=1980+3.7*T0*u[5] 
> EE=(9.0*N0-2.0*u[1]-u[2])/(9.0*N0-6.0*u[1]-3.0*u[2])
> dEE=18.0*N0/(9.0*N0-6.0*u[1]-3.0*u[2])^2
> lambda=lambda_T0*exp(Ea/(8.314*T0)*((1.0/u[5])-1.0)) 
> Cf=0.205*((ro_a*u[4]*u0*u[3]*R0)/eta_a)^(-0.61) 
> C2=2.0*Cf*ro_a*L*u0^2/(G0*R0)
> 
> M = [-u[4] 0 0 -u[1] 0
> 0 -u[4] 0 +2.0*u[2] 0
> 0 0 -u[4] -u[3]/2.0 0
> -2.0*dEE*(u[2]-u[1])+EE -(dEE*(u[2]-u[1])+EE) -2.0*EE*(u[2]-u[1])/u[3] +C1*u[4] 0
> 
> 0 0 0 -V3*EE*(u[2]-u[1])/(Cp*u[4]) +1.0]
> 
> end
> 
> u₀ = [1.0/Eeq, czz0, 1.0, 1.0, 1.0] 
> tspan = (0.0,1.0)
>     
> Cp0=1980+3.7*T0*u₀[5] #969.9 
> EE0=(9.0*N0-2.0*u₀[1]-u₀[2])/(9.0*N0-6.0*u₀[1]-3.0*u₀[2]) 
> dEE0=18.0*N0/(9.0*N0-6.0*u₀[1]-3.0*u₀[2])^2 
> lambda0=lambda_T0*exp(Ea/(8.314*T0)*((1.0/u₀[5])-1.0)) 
> Cf0=0.205*((ro_a*u₀[4]*u0*u₀[3]*R0)/eta_a)^(-0.61) 
> C20=2*Cf0*ro_a*L*u0^2/(G0*R0)      
> 
> M0 = [-u₀[4] 0 0 -u₀[1] 0 
> 0 -u₀[4] 0 +2.0*u₀[2] 0 
> 0 0 -u₀[4] -u₀[3]/2 0 
> -2.0*dEE0*(u₀[2]-u₀[1])+EE0 -(dEE0*(u₀[2]-u₀[1])+EE0) -2.0*EE0*(u₀[2]-u₀[1])/u₀[3] +C1*u₀[4] 0 
> 0 0 0 -V3*EE0*(u₀[2]-u₀[1])/(Cp0*u₀[4]) +1]
> 
> M = DiffEqArrayOperator(M0,update_func=update_func)
> 
> f = ODEFunction(filament,mass_matrix=M)
> prob_mm = ODEProblem(f,u₀,tspan)
> 
> sol = DifferentialEquations.solve(prob_mm, RadauIIA5())
> 
> display(plot(sol, tspan=(0.0, 1.0), layout=(5,1)))
> 
> ```

Besides, there is a message

> Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(filament),DiffEqArrayOperator{Float64,Array{Float64,2},typeof(update\_func)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64}}(1.0,0.0,0.0,0.0,1.0,0.0)Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(filament),DiffEqArrayOperator{Float64,Array{Float64,2},typeof(update\_func)},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64}}(1.0,0.0,0.0,1.0,0.0,0.0)

Is something wrong here?

---

<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:** [July 17, 2020, 7:02am UTC](https://discourse.julialang.org/t/dependent-mass-matrix-with-radauiia5/43189/2 "2020-07-17T07:02:17Z")

</div>

Dependent mass matrices like that are pretty fundamentally incompatible with automatic differentiation (if you want to do it in a way that’s fast), so after experimenting a lot with this, we are going to be dropping the experimental support for dependent mass matrices. That said, any dependent mass matrix problem can be rewritten as a constant mass matrix problem, which is a numerically much better approach (and something we hope to automate through ModelingToolkit fairly soon).

---

<div class="post-metadata">

**Author:** ![Jorge\_Silva](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jorge_silva/32/14976_2.png) [@Jorge\_Silva](https://discourse.julialang.org/u/Jorge_Silva)\
**Post date:** [July 17, 2020, 2:56pm UTC](https://discourse.julialang.org/t/dependent-mass-matrix-with-radauiia5/43189/3 "2020-07-17T14:56:42Z")

</div>

Thanks Chris.  
If I understood correctly, converting the system to the form du=M-1(t,u)f(t,u) will give better reuslts.  
In addition I ve rewritten the system as f(t,u,du)=0 and runned Sundials IDA; the result is good for some conditions. However, depending on the value of parameters the solution can diverge very quickly.

---

<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:** [July 17, 2020, 3:33pm UTC](https://discourse.julialang.org/t/dependent-mass-matrix-with-radauiia5/43189/4 "2020-07-17T15:33:10Z")

</div>

> [@Jorge\_Silva](#):
>
> converting the system to the form du=M-1(t,u)f(t,u) will give better reuslts

No, don’t invert it (that won’t be possible with DAEs anyways). Just write the version that’s the equivalent equation with constant mass matrices through substitution.

---

<div class="post-metadata">

**Author:** ![Jorge\_Silva](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jorge_silva/32/14976_2.png) [@Jorge\_Silva](https://discourse.julialang.org/u/Jorge_Silva)\
**Post date:** [July 17, 2020, 5:28pm UTC](https://discourse.julialang.org/t/dependent-mass-matrix-with-radauiia5/43189/5 "2020-07-17T17:28:48Z")

</div>

Got it! I going to try that way.  
Thanks