Error while using DifferentialEquations.jl stiff methods

#1

Hi all,

The DifferentialEquations’ suggested integration methods when trying to solve this stiff problem all yield the same error. The native Julia ones, that is. Sundial’s CVODE_BDF as well as lsoda are able to solve this with no problems.

My code:

using DifferentialEquations
using Plots
gr()
using LSODA
using Sundials

x0 = zeros(Float64, 35)
x0[1] = 5
x0[2] = 1
x0[3] = 5
x0[4] = 5
x0[5] = 1
x0[6] = 1
x0[7] = 1
x0[8] = 1
x0[9] = 1
x0[10] = 1
x0[11] = 1
x0[12] = 1
x0[13] = 1
x0[14] = 1
x0[15] = 1
x0[16] = 1
x0[17] = 1
x0[18] = 1
x0[19] = 1
x0[20] = 1
x0[21] = 1
x0[22] = 1
x0[23] = 1
x0[24] = 1
x0[25] = 1
x0[26] = 1
x0[27] = 1
x0[28] = 1
x0[29] = 1
x0[30] = 1
x0[31] = 1
x0[32] = 1
x0[33] = 1
x0[34] = 1
x0[35] = 6

p = zeros(Float64, 117)
p[1] = 1000
p[2] = 1000
p[3] = 1000
p[4] = 1000
p[5] = 1000
p[6] = 1000
p[7] = 1000
p[8] = 1000
p[9] = 1
p[10] = 1
p[11] = 0.7
p[12] = 0.7
p[13] = 0.7
p[14] = 0.7
p[15] = 0.2
p[16] = 0.2
p[17] = 1
p[18] = 1
p[19] = 0.7
p[20] = 0.7
p[21] = 0.05
p[22] = 1
p[23] = 0.5
p[24] = 0.7
p[25] = 0.5
p[26] = 2
p[27] = 0.7
p[28] = 0.7
p[29] = 2
p[30] = 0.2
p[31] = 0.2
p[32] = 1
p[33] = 0.7
p[34] = 0.2
p[35] = 0.21
p[36] = 1
p[37] = 0.7
p[38] = 1
p[39] = 0.2
p[40] = 0.2
p[41] = 0.05
p[42] = 0.01
p[43] = 0.7
p[44] = 1
p[45] = 0.5
p[46] = 1
p[47] = 0.04
p[48] = 0.01
p[49] = 0.01
p[50] = 0.01
p[51] = 1
p[52] = 1
p[53] = 0.5
p[54] = 1
p[55] = 1
p[56] = 0.5
p[57] = 0.1
p[58] = 0.1
p[59] = 0.5
p[60] = 1
p[61] = 1
p[62] = 1
p[63] = 0.3
p[64] = 1
p[65] = 1
p[66] = 0.7
p[67] = 0.4
p[68] = 0.7
p[69] = 1
p[70] = 0.7
p[71] = 0.7
p[72] = 2
p[73] = 0.2
p[74] = 2
p[75] = 1
p[76] = 0.7
p[77] = 0.7
p[78] = 0.2
p[79] = 0.7
p[80] = 0.5
p[81] = 0.21
p[82] = 0.2
p[83] = 0.2
p[84] = 0.01
p[85] = 2
p[86] = 2
p[87] = 2
p[88] = 2
p[89] = 2
p[90] = 2
p[91] = 2
p[92] = 2
p[93] = 1
p[94] = 1
p[95] = 1
p[96] = 1
p[97] = 1
p[98] = 0.5
p[99] = 0.7
p[100] = 0.7
p[101] = 2
p[102] = 0.2
p[103] = 2
p[104] = 1
p[105] = 1
p[106] = 0.7
p[107] = 1
p[108] = 0.7
p[109] = 1
p[110] = 0.7
p[111] = 1
p[112] = 0.7
p[113] = 2
p[114] = 2
p[115] = 0.7
p[116] = 0.2
p[117] = 0.2

function fun_b4(dx_dt, x, par, t)
    p = zeros(347)
    p[118] = 1
    p[119] = 1
    p[120] = 1
    p[121] = 1
    p[122] = 1
    p[123] = 1
    p[124] = 1
    p[125] = 1
    p[126] = 1
    p[127] = 1
    p[128] = 1
    p[129] = 1
    p[130] = 1
    p[131] = 1
    p[132] = 1
    p[133] = 1
    p[134] = 1
    p[135] = 1
    p[136] = 1
    p[137] = 1
    p[138] = 1
    p[139] = 1
    p[140] = 1
    p[141] = 1
    p[142] = 1
    p[143] = 1
    p[144] = 1
    p[145] = 1
    p[146] = 1
    p[147] = 1
    p[148] = 1
    p[149] = 1
    p[150] = 0
    p[151] = 141.471
    p[152] = 0.9
    p[153] = 0.1
    p[154] = 83402
    p[155] = 127877
    p[156] = 603.414
    p[157] = 603.414
    p[158] = 5.02845
    p[159] = 21722.9
    p[160] = 603.414
    p[161] = 361200
    p[162] = 673524
    p[163] = 1.62679e+06
    p[164] = 164390
    p[165] = 36512.9
    p[166] = 166804
    p[167] = 39530
    p[168] = 1810.24
    p[169] = 127877
    p[170] = 166201
    p[171] = 1.94454e+06
    p[172] = 5.02845
    p[173] = 317755
    p[174] = 0
    p[175] = 532531
    p[176] = 38926.6
    p[177] = 40133.4
    p[178] = 689982
    p[179] = 4.45873e+06
    p[180] = 2413.66
    p[181] = 83402
    p[182] = 603.414
    p[183] = 603.414
    p[184] = 2413.66
    p[185] = 533135
    p[186] = 100000
    p[187] = 1
    p[188] = 1000
    p[189] = 1000
    p[190] = 30
    p[191] = 1000
    p[192] = 1000
    p[193] = 1000
    p[194] = 1000
    p[195] = 1000
    p[196] = 1000
    p[197] = 1000
    p[198] = 3000
    p[199] = 1000
    p[200] = 1000
    p[201] = 1000
    p[202] = 100
    p[203] = 100
    p[204] = 100
    p[205] = 100
    p[206] = 1000
    p[207] = 1000
    p[208] = 1000
    p[209] = 1000
    p[210] = 1000
    p[211] = 1000
    p[212] = 1000
    p[213] = 1000
    p[214] = 1000
    p[215] = 3000
    p[216] = 100000
    p[217] = 1000
    p[218] = 1000
    p[219] = 1000
    p[220] = 1000
    p[221] = 1000
    p[222] = 1000
    p[223] = 1000
    p[224] = 1000
    p[225] = 1000
    p[226] = 1000
    p[227] = 1000
    p[228] = 0
    p[229] = 1
    p[230] = 0
    p[231] = 0
    p[232] = 1
    p[233] = 1
    p[234] = 1
    p[235] = -0.7
    p[236] = -0.7
    p[237] = 0.7
    p[238] = 0.7
    p[239] = -0.2
    p[240] = -0.2
    p[241] = 1
    p[242] = 1
    p[243] = -0.7
    p[244] = -0.7
    p[245] = 0.05
    p[246] = 1
    p[247] = 0.5
    p[248] = -0.7
    p[249] = -0.5
    p[250] = 2
    p[251] = 0.7
    p[252] = 0.7
    p[253] = -2
    p[254] = -0.2
    p[255] = -0.2
    p[256] = 1
    p[257] = 0.7
    p[258] = -0.2
    p[259] = -0.21
    p[260] = 0
    p[261] = 1
    p[262] = 0.7
    p[263] = 1
    p[264] = -0.2
    p[265] = -0.2
    p[266] = 0.05
    p[267] = 0
    p[268] = -0.01
    p[269] = 0.7
    p[270] = 1
    p[271] = -0.5
    p[272] = -1
    p[273] = 0.04
    p[274] = -0.01
    p[275] = -0.01
    p[276] = -0.01
    p[277] = 1
    p[278] = 1
    p[279] = 0.5
    p[280] = -1
    p[281] = -1
    p[282] = -0.5
    p[283] = 0.1
    p[284] = -0.1
    p[285] = 0.5
    p[286] = 1
    p[287] = 1
    p[288] = -1
    p[289] = -0.3
    p[290] = -1
    p[291] = 1
    p[292] = 0.7
    p[293] = -0.4
    p[294] = -0.7
    p[295] = 1
    p[296] = -0.7
    p[297] = 0.7
    p[298] = 2
    p[299] = -0.2
    p[300] = -2
    p[301] = 1
    p[302] = 0.7
    p[303] = 0.7
    p[304] = 0.2
    p[305] = 0.7
    p[306] = -0.5
    p[307] = -0.21
    p[308] = -0.2
    p[309] = -0.2
    p[310] = 0
    p[311] = 0
    p[312] = -0.01
    p[313] = 0.7
    p[314] = -0.2
    p[315] = 2
    p[316] = 2
    p[317] = -2
    p[318] = -2
    p[319] = 2
    p[320] = 2
    p[321] = -2
    p[322] = -2
    p[323] = 1
    p[324] = 1
    p[325] = -1
    p[326] = -1
    p[327] = 1
    p[328] = -0.5
    p[329] = 0.7
    p[330] = 0.7
    p[331] = 2
    p[332] = -0.2
    p[333] = -2
    p[334] = 1
    p[335] = 1
    p[336] = -0.7
    p[337] = 1
    p[338] = -0.7
    p[339] = 1
    p[340] = -0.7
    p[341] = 1
    p[342] = -0.7
    p[343] = 2
    p[344] = -2
    p[345] = 0.7
    p[346] = -0.2
    p[347] = -0.2

    for i in 1:length(par)
        p[i] = par[i]
    end
    #endif /* FIXED */

    #ifdef ASSIGNMENT
    y = zeros(102)
    y[1] = log(x[1])
    y[2] = log(x[2])
    y[3] = log(x[3])
    y[4] = log(x[4])
    y[5] = log(x[5])
    y[6] = log(x[6])
    y[7] = log(x[7])
    y[8] = log(x[8])
    y[9] = log(x[9])
    y[10] = log(x[10])
    y[11] = log(x[11])
    y[12] = log(x[12])
    y[13] = log(x[13])
    y[14] = log(x[14])
    y[15] = log(x[15])
    y[16] = log(x[16])
    y[17] = log(x[17])
    y[18] = log(x[18])
    y[19] = log(x[19])
    y[20] = log(x[20])
    y[21] = log(x[21])
    y[22] = log(x[22])
    y[23] = log(x[23])
    y[24] = log(x[24])
    y[25] = log(x[25])
    y[26] = log(x[26])
    y[27] = log(x[27])
    y[28] = log(x[28])
    y[29] = log(x[29])
    y[30] = log(x[30])
    y[31] = log(x[31])
    y[32] = log(x[32])
    y[33] = log(x[33])
    y[34] = log(x[34])
    y[35] = log(x[35])
    y[36] = x[1]*p[186]
    y[37] = x[2]*p[187]
    y[38] = x[3]*p[188]
    y[39] = x[4]*p[189]
    y[40] = x[5]*p[190]
    y[41] = x[6]*p[191]
    y[42] = x[7]*p[192]
    y[43] = x[8]*p[193]
    y[44] = x[9]*p[194]
    y[45] = x[10]*p[195]
    y[46] = x[11]*p[196]
    y[47] = x[12]*p[197]
    y[48] = x[13]*p[198]
    y[49] = x[14]*p[199]
    y[50] = x[15]*p[200]
    y[51] = x[16]*p[201]
    y[52] = x[17]*p[202]
    y[53] = x[18]*p[203]
    y[54] = x[19]*p[204]
    y[55] = x[20]*p[205]
    y[56] = x[21]*p[206]
    y[57] = x[22]*p[207]
    y[58] = x[23]*p[208]
    y[59] = x[24]*p[209]
    y[60] = x[25]*p[210]
    y[61] = x[26]*p[211]
    y[62] = x[27]*p[212]
    y[63] = x[28]*p[213]
    y[64] = x[29]*p[214]
    y[65] = x[30]*p[215]
    y[66] = x[31]*p[216]
    y[67] = x[32]*p[217]
    y[68] = x[33]*p[218]
    y[69] = x[34]*p[219]
    if t>5 && t<10
        y[70] = p[150]
    else
        y[70] = 0.00000000000000000
    end
    y[71] = 0.00000000000000000
    y[72] = 0.00000000000000000
    y[73] = 0.00000000000000000
    y[74] = 0.00000000000000000
    y[75] = 0.00000000000000000
    y[76] = p[123]*p[159]*(1.00000000000000000+p[233]*y[8]+p[234]*y[9]+p[235]*y[6]+p[236]*y[7])
    y[77] = p[124]*p[160]*(1.00000000000000000+p[237]*y[12]+p[238]*y[13]+p[239]*y[10]+p[240]*y[11])
    y[78] = p[125]*p[161]*(1.00000000000000000+p[241]*y[15]+p[242]*y[6]+p[243]*y[8]+p[244]*y[14]+p[245]*y[16])
    y[79] = p[126]*p[162]*(1.00000000000000000+p[246]*y[16]+p[247]*y[7]+p[248]*y[15]+p[249]*y[9])
    y[80] = p[127]*p[163]*(1.00000000000000000+p[250]*y[9]+p[251]*y[19]+p[252]*y[20]+p[253]*y[7]+p[254]*y[17]+p[255]*y[18])
    y[81] = p[128]*p[164]*(1.00000000000000000+p[256]*y[22]+p[257]*y[11]+p[258]*y[21]+p[259]*y[13]+p[260]*y[13])
    y[82] = p[129]*p[165]*(1.00000000000000000+p[261]*y[15]+p[262]*y[7]+p[263]*y[21]+p[264]*y[8]+p[265]*y[9]+p[266]*y[32]+p[267]*y[9]+p[268]*y[30])
    y[83] = p[130]*p[166]*(1.00000000000000000+p[269]*y[25]+p[270]*y[26]+p[271]*y[23]+p[272]*y[24]+p[273]*y[11]+p[274]*y[11]+p[275]*y[13]+p[276]*y[22])
    y[84] = p[131]*p[167]*(1.00000000000000000+p[277]*y[28]+p[278]*y[29]+p[279]*y[23]+p[280]*y[12]+p[281]*y[27]+p[282]*y[25]+p[283]*y[27]+p[284]*y[10])
    y[85] = p[132]*p[168]*(1.00000000000000000+p[285]*y[32]+p[286]*y[22]+p[287]*y[16]+p[288]*y[15]+p[289]*y[30]+p[290]*y[27])
    y[86] = p[133]*p[169]*(1.00000000000000000+p[291]*y[21]+p[292]*y[23]+p[293]*y[2]+p[294]*y[25])
    y[87] = p[134]*p[170]*(1.00000000000000000+p[295]*y[24]+p[296]*y[22])
    y[88] = p[135]*p[171]*(1.00000000000000000+p[297]*y[20]+p[298]*y[17]+p[299]*y[18]+p[300]*y[19])
    y[89] = p[136]*(p[172]/(p[209]/p[220]*(p[210]/p[221])*(p[197]/p[222])*(p[218]/p[223])*(p[219]/p[224])*(p[214]/p[225])*(p[195]/p[226])*(p[198]/p[227])/((1.00000000000000000+p[209]/p[220])*(1.00000000000000000+p[210]/p[221])*(1.00000000000000000+p[197]/p[222])*(1.00000000000000000+p[218]/p[223])*(1.00000000000000000+p[219]/p[224])*(1.00000000000000000+p[214]/p[225])*(1.00000000000000000+p[195]/p[226])*(1.00000000000000000+p[198]/p[227]))))*(x[24]*p[209]/p[220]*(x[25]*p[210]/p[221])*(x[12]*p[197]/p[222])*(x[33]*p[218]/p[223])*(x[34]*p[219]/p[224])*(x[29]*p[214]/p[225])*(x[10]*p[195]/p[226])*(x[13]*p[198]/p[227])/((1.00000000000000000+x[24]*p[209]/p[220])*(1.00000000000000000+x[25]*p[210]/p[221])*(1.00000000000000000+x[12]*p[197]/p[222])*(1.00000000000000000+x[33]*p[218]/p[223])*(1.00000000000000000+x[34]*p[219]/p[224])*(1.00000000000000000+x[29]*p[214]/p[225])*(1.00000000000000000+x[10]*p[195]/p[226])*(1.00000000000000000+x[13]*p[198]/p[227])))
    y[90] = p[137]*p[173]*(1.00000000000000000+p[301]*y[8]+p[302]*y[7]+p[303]*y[19]+p[304]*y[31]+p[305]*y[32]+p[306]*y[16]+p[307]*y[9]+p[308]*y[17]+p[309]*y[30]+p[310]*y[31]+p[311]*y[9]+p[312]*y[15])
    y[91] = p[138]*p[174]*(1.00000000000000000+p[313]*y[18]+p[314]*y[20])
    y[92] = p[139]*p[175]*(1.00000000000000000+p[315]*y[11]+p[316]*y[30]+p[317]*y[32]+p[318]*y[13])
    y[93] = p[140]*p[176]*(1.00000000000000000+p[319]*y[8]+p[320]*y[27]+p[321]*y[28]+p[322]*y[16])
    y[94] = p[141]*p[177]*(1.00000000000000000+p[323]*y[14]+p[324]*y[12]+p[325]*y[29]+p[326]*y[6])
    y[95] = p[142]*p[178]*(1.00000000000000000+p[327]*y[13]+p[328]*y[11])
    y[96] = p[143]*p[179]*(1.00000000000000000+p[329]*y[32]+p[330]*y[31]+p[331]*y[18]+p[332]*y[30]+p[333]*y[20])
    y[97] = p[144]*p[180]*(1.00000000000000000+p[334]*y[27]+p[335]*y[31]+p[336]*y[16])
    y[98] = p[145]*p[181]*(1.00000000000000000+p[337]*y[1]+p[338]*y[26])
    y[99] = p[146]*p[182]*(1.00000000000000000+p[339]*y[3]+p[340]*y[33])
    y[100] = p[147]*p[183]*(1.00000000000000000+p[341]*y[4]+p[342]*y[34])
    y[101] = p[148]*p[184]*(1.00000000000000000+p[343]*y[6]+p[344]*y[12])
    y[102] = p[149]*p[185]*(1.00000000000000000+p[345]*y[18]+p[346]*y[31]+p[347]*y[20])
    #endif /* ASSIGNMENT */

    #ifdef ODEs
    dx_dt[1] = (y[71]-y[98]*1.00000000000000000/p[151]-x[1]*p[186]/x[35]*y[70])/p[186]
    dx_dt[2] = ((-y[72])+y[86]*1.00000000000000000/p[151]-x[2]*p[187]/x[35]*y[70])/p[187]
    dx_dt[3] = (y[73]-y[99]*1.00000000000000000/p[151]-x[3]*p[188]/x[35]*y[70])/p[188]
    dx_dt[4] = (y[74]-y[100]*1.00000000000000000/p[151]-x[4]*p[189]/x[35]*y[70])/p[189]
    dx_dt[5] = ((-y[75])+y[89]*1.00000000000000000/p[151]-x[5]*p[190]/x[35]*y[70])/p[190]
    dx_dt[6] = (y[76]-y[78]+y[94]*p[152]/p[153]-y[101]*p[152]/p[153])/p[191]
    dx_dt[7] = (y[76]-y[79]+y[80]-2.00000000000000000*y[82]*p[152]/p[153]-y[90])/p[192]
    dx_dt[8] = ((-y[76])+y[78]+y[82]*p[152]/p[153]-y[90]-y[93]*p[152]/p[153])/p[193]
    dx_dt[9] = ((-y[76])+y[79]-y[80]+2.00000000000000000*y[82]*p[152]/p[153]+y[90])/p[194]
    dx_dt[10] = (y[77]-120.00000000000000000*y[89])/p[195]
    dx_dt[11] = (y[77]-y[81]+1260.00000000000000000*y[89]-y[92]+y[95])/p[196]
    dx_dt[12] = ((-y[77])+y[84]-240.00000000000000000*y[89]-y[94]+y[101])/p[197]
    dx_dt[13] = ((-y[77])+y[81]-1260.00000000000000000*y[89]+y[92]-y[95])/p[198]
    dx_dt[14] = (y[78]-y[94]*p[152]/p[153])/p[199]
    dx_dt[15] = ((-y[78])+y[79]-y[82]*p[152]/p[153]+y[85]*p[152]/p[153])/p[200]
    dx_dt[16] = ((-y[79])-y[85]*p[152]/p[153]+y[90]+y[93]*p[152]/p[153]+y[97]*p[152]/p[153])/p[201]
    dx_dt[17] = (2.00000000000000000*y[80]-2.00000000000000000*y[88]+2.00000000000000000*y[90])/p[202]
    dx_dt[18] = (4.00000000000000000*y[80]+6.00000000000000000*y[88]-y[91]-3.00000000000000000*y[96]-y[102]*p[152]/p[153])/p[203]
    dx_dt[19] = ((-2.00000000000000000)*y[80]+2.00000000000000000*y[88]-2.00000000000000000*y[90])/p[204]
    dx_dt[20] = ((-4.00000000000000000)*y[80]-6.00000000000000000*y[88]+y[91]+3.00000000000000000*y[96]+y[102]*p[152]/p[153])/p[205]
    dx_dt[21] = (y[81]-y[82]-y[86])/p[206]
    dx_dt[22] = ((-y[81])-y[85]+y[87])/p[207]
    dx_dt[23] = (y[83]-y[84]-y[86]+120.00000000000000000*y[89])/p[208]
    dx_dt[24] = (y[83]-y[87]-120.00000000000000000*y[89])/p[209]
    dx_dt[25] = ((-y[83])+y[84]+y[86]-120.00000000000000000*y[89])/p[210]
    dx_dt[26] = ((-0.50000000000000000)*y[83]+y[98])/p[211]
    dx_dt[27] = (y[84]+y[85]-y[93]-y[97])/p[212]
    dx_dt[28] = ((-y[84])+120.00000000000000000*y[89]+y[93])/p[213]
    dx_dt[29] = ((-y[84])-120.00000000000000000*y[89]+y[94])/p[214]
    dx_dt[30] = (y[85]*p[152]/p[153]+y[90]-y[92]*p[152]/p[153]+y[96])/p[215]
    dx_dt[31] = ((-y[90])-y[96]-y[97]*p[152]/p[153]+y[102]*p[152]/p[153])/p[216]
    dx_dt[32] = ((-y[85])*p[152]/p[153]-y[90]+y[92]*p[152]/p[153]-y[96])/p[217]
    dx_dt[33] = ((-120.00000000000000000)*y[89]+y[99])/p[218]
    dx_dt[34] = ((-120.00000000000000000)*y[89]+y[100])/p[219]
    dx_dt[35] = y[70]
end



tspan = (1.0, 20.0)
solvers_list = [lsoda(), CVODE_BDF(), Rosenbrock23(), TRBDF2(), ABDF2(), Rodas5(), Rodas4P(), Kvaerno5(), KenCarp4(), RadauIIA5()]

for solver in solvers_list
    try
        prob = ODEProblem(fun_b4, x0, tspan, p)
        sol  = solve(prob, solver)
        plot_sol = plot(sol, title="$solver"[1:match(r"(\{)|(\()", "$solver").offset-1])
        display(plot_sol)
    catch e
        @warn("\nException:\n$e\nFor solver:\n$solver\n")
    end
end

#=>
# One plot for each chemical element
sol_u = reduce(hcat, sol.u)
for i in 1:size(sol_u)[1]
    p = plot(sol.t, sol_u[i,:])
    display(p)
end
<=#

And the errors:

┌ Warning: 
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ Rosenbrock23{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), Val{:central})
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
┌ Warning: 
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ TRBDF2{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},NLNewton{Nothing,Nothing},DataType,Float64,:Predictive}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), NLNewton{Nothing,Nothing}(nothing, nothing, 10), Val{:central}, true, :linear, 0.001)
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
┌ Warning: 
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ ABDF2{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},NLNewton{Nothing,Nothing},DataType,Nothing,Nothing,Float64,:Predictive}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), NLNewton{Nothing,Nothing}(nothing, nothing, 10), Val{:central}, nothing, nothing, true, :linear, 0.001)
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
┌ Warning: 
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ Rodas5{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), Val{:central})
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
┌ Warning: 
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ Rodas4P{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), Val{:central})
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
┌ Warning: 
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ Kvaerno5{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},NLNewton{Nothing,Nothing},DataType,Float64,:Predictive}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), NLNewton{Nothing,Nothing}(nothing, nothing, 10), Val{:central}, true, :linear, 0.001)
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
┌ Warning: 
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ KenCarp4{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},NLNewton{Nothing,Nothing},DataType,Float64,:Predictive}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), NLNewton{Nothing,Nothing}(nothing, nothing, 10), Val{:central}, true, :linear, 0.001)
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
┌ Warning: 
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ RadauIIA5{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType,Float64,Nothing,:Predictive}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), Val{:central}, true, :dense, 0.001, nothing, nothing, 10, 1)
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567

Not sure if this is an actual issue or if I’m doing something wrong here.

As a side note: in case anyone knows a better way for plotting each chemical element, as mentioned in the last comment of my code, please let me know. :slight_smile:


– Edit:
I just ran this code again and this time lsoda got lost as well.

┌ Warning: 
│ Exception:
│ DomainError(-0.14279345102360447, "log will only return a complex result if called with a complex argument. Try log(Complex(x)).")
│ For solver:
│ lsoda()
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:561

So the error here is clearer for me, it’s having trouble with the log function.
But it seems like most of the time lsoda will get it right, with this error popping up just a couple of times after some 20+ runs. CVODE_BDF doesn’t seem to have a problem with log at all.

#2

Your issue is that you’re not handling autodiff’d values correctly.

That assumes that all of your values will be Float64, which isn’t true during autodiff. So either turn autodiff off like Rosenbrock23(autodiff=false), or make sure your types are correct, like:

y = zeros(promote_type(eltype(x),typeof(t)),102)
1 Like
#3

I see!
Thanks @ChrisRackauckas, I’ll be sure to keep my typing in mind from now on.

#4

No worries, getting used to autodiff is quite a step, but it has a lot of performance and numerical stability benefits so it’s a good thing to try. But the stiff solvers that use it have an argument that is autodiff=false to remove it and resort back to numerical differentiation if necessary.