From 5953e551837a56c60ac04308c5c669cf368ef1a0 Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Thu, 20 Mar 2014 16:23:32 +0100 Subject: [PATCH 01/20] Initial import of dop853. --- src/ODE.jl | 3 +- src/dop853.jl | 418 ++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 420 insertions(+), 1 deletion(-) create mode 100644 src/dop853.jl diff --git a/src/ODE.jl b/src/ODE.jl index 116a0fcbe..d6e0d9e56 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -5,7 +5,7 @@ module ODE using Polynomial ## minimal function export list -export ode23, ode4, ode45, ode4s, ode4ms +export ode23, ode4, ode45, ode4s, ode4ms, dop853 ## complete function export list #export ode23, ode4, @@ -13,6 +13,7 @@ export ode23, ode4, ode45, ode4s, ode4ms # oderosenbrock, ode4s, ode4s_kr, ode4s_s, # ode4ms, ode_ms +include("dop853.jl") #ODE23 Solve non-stiff differential equations. # diff --git a/src/dop853.jl b/src/dop853.jl new file mode 100644 index 000000000..9309dcdf4 --- /dev/null +++ b/src/dop853.jl @@ -0,0 +1,418 @@ +function dop853(f::Function, x::Float64, yin::Vector{Float64}, xend::Float64, params::Any; rtol::Vector{Float64}=[1e-6], atol::Vector{Float64}=[sqrt(eps())], + uround::Float64=2.3e-16, safe::Float64=0.9, fac1::Float64=0.333, fac2::Float64=6.0, + beta::Float64=0.0, maxstepsize::Float64=xend-x, initialstep::Float64=0.0, + maxsteps::Int64=100000, printmessages::Bool=false, nstiff::Int64=1000, + iout::Int64=0, solout::Function=s(x...)=return, dense::Vector{Int64}=[1:length(yin)]) + + c14 = 0.1e+00 + c15 = 0.2e+00 + c16 = 0.777777777777777777777777777778e+00 + d41 = -0.84289382761090128651353491142e+01 + d46 = 0.56671495351937776962531783590e+00 + d47 = -0.30689499459498916912797304727e+01 + d48 = 0.23846676565120698287728149680e+01 + d49 = 0.21170345824450282767155149946e+01 + d410 = -0.87139158377797299206789907490e+00 + d411 = 0.22404374302607882758541771650e+01 + d412 = 0.63157877876946881815570249290e+00 + d413 = -0.88990336451333310820698117400e-01 + d414 = 0.18148505520854727256656404962e+02 + d415 = -0.91946323924783554000451984436e+01 + d416 = -0.44360363875948939664310572000e+01 + d51 = 0.10427508642579134603413151009e+02 + d56 = 0.24228349177525818288430175319e+03 + d57 = 0.16520045171727028198505394887e+03 + d58 = -0.37454675472269020279518312152e+03 + d59 = -0.22113666853125306036270938578e+02 + d510 = 0.77334326684722638389603898808e+01 + d511 = -0.30674084731089398182061213626e+02 + d512 = -0.93321305264302278729567221706e+01 + d513 = 0.15697238121770843886131091075e+02 + d514 = -0.31139403219565177677282850411e+02 + d515 = -0.93529243588444783865713862664e+01 + d516 = 0.35816841486394083752465898540e+02 + d61 = 0.19985053242002433820987653617e+02 + d66 = -0.38703730874935176555105901742e+03 + d67 = -0.18917813819516756882830838328e+03 + d68 = 0.52780815920542364900561016686e+03 + d69 = -0.11573902539959630126141871134e+02 + d610 = 0.68812326946963000169666922661e+01 + d611 = -0.10006050966910838403183860980e+01 + d612 = 0.77771377980534432092869265740e+00 + d613 = -0.27782057523535084065932004339e+01 + d614 = -0.60196695231264120758267380846e+02 + d615 = 0.84320405506677161018159903784e+02 + d616 = 0.11992291136182789328035130030e+02 + d71 = -0.25693933462703749003312586129e+02 + d76 = -0.15418974869023643374053993627e+03 + d77 = -0.23152937917604549567536039109e+03 + d78 = 0.35763911791061412378285349910e+03 + d79 = 0.93405324183624310003907691704e+02 + d710 = -0.37458323136451633156875139351e+02 + d711 = 0.10409964950896230045147246184e+03 + d712 = 0.29840293426660503123344363579e+02 + d713 = -0.43533456590011143754432175058e+02 + d714 = 0.96324553959188282948394950600e+02 + d715 = -0.39177261675615439165231486172e+02 + d716 = -0.14972683625798562581422125276e+03 + + irtrn = 0 + + n = length(yin) + nrd = length(dense) + nfcn = 0 + nstep = 0 + naccpt = 0 + nrejct = 0 + k1 = zeros(n) + k2 = zeros(n) + k3 = zeros(n) + k4 = zeros(n) + k5 = zeros(n) + k6 = zeros(n) + k7 = zeros(n) + k8 = zeros(n) + k9 = zeros(n) + k10 = zeros(n) + y = zeros(n) + copy!(y, yin) + facold = 1e-4 + expo1 = 1.0/8.0 - beta*0.2 + facc1 = 1.0/fac1 + facc2 = 1.0/fac2 + posneg = sign(xend-x) + atol = length(atol) == 1 ? fill(atol[1], n) : atol + rtol = length(rtol) == 1 ? fill(rtol[1], n) : rtol + last = false + hlamb = 0.0 + iasti = 0 + f(k1, x, y, params) + nfcn += 1 + hmax = abs(maxstepsize) + nmax = maxsteps + h = initialstep + iord = 8 + if h == 0.0 + h = hinit(n, f, x, y, xend, posneg, k1, k2, k3, iord, hmax, atol, rtol, params) + nfcn += 1 + end + reject = false + xold = x + if iout != 0 + #hout = 1.0 + irtrn = solout(naccpt+1, xold, x, y, con, icomp, params) + if irtrn < 0 + return y + end + end + + while true + #if nstep == 4 + # println(y) + # exit() + #end + if nstep > nmax + error("Exit at x=$x. More than nmax=$nmax steps needed.") + end + if 0.1*abs(h) <= abs(x)*uround + error("Exit at x=$x. Step size too small, h=$h.") + end + if (x+1.01*h-xend)*posneg > 0.0 + h = xend-x + last = true + end + nstep += 1 + if irtrn >= 2 + f(k1, x, y, params) + end + + y1, err = dopcore(nstep, f, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, atol, rtol, params) + xph = x+h + nfcn += 11 + + fac11 = err^expo1 + fac = fac11/facold^beta + fac = max(facc2, min(facc1, fac/safe)) + hnew = h/fac + if err <= 1.0 + facold = max(err, 1e-4) + naccpt += 1 + f(k4, xph, k5, params) + nfcn += 1 + # Stiffness detection + if mod(naccpt, nstiff) == 0 || iasti > 0 + nonsti = 0 + stnum = 0.0 + stden = 0.0 + for i = 1:n + stnum += (k4[i] - k3[i])^2 + stden += (k5[i] - y1[i])^2 + end + if stden > 0.0 + hlamb = abs(h)*sqrt(stnum/stden) + end + if hlamb > 6.1 + nonsti = 0 + iasti += 1 + if iasti == 15 + if printmessages + println("The problem seems to become stiff at x=$x") + end + end + else + nonsti += 1 + if nonsti == 6 + iasti = 0 + end + end + end + # Final preparation for dense output + event = iout == 3 && xout <= xph + if iout == 2 || event + for (i,j) in enumerate(dense) + end + end + copy!(k1, k4) + copy!(y, k5) + xold = x + x = xph + # if + # solout + # end + # Normal exit + if last + h = hnew + idid = 1 + return y + end + if abs(hnew) > hmax + hnew = posneg*hmax + end + if reject + hnew = posneg*min(abs(hnew),abs(h)) + end + reject = false + else + hnew = h/min(facc1,fac11/safe) + reject = true + if naccpt >= 1 + nrejct += 1 + end + last = false + end + h = hnew + end +end + +function hinit(n::Int64, f::Function, x::Float64, y::Vector{Float64}, xend::Float64, posneg::Float64, f0::AbstractArray{Float64,1}, f1::AbstractArray{Float64,1}, y0::AbstractArray{Float64,1}, iord::Int64, hmax::Float64, atol::Vector{Float64}, rtol::Vector{Float64}, params::Any) + dnf = 0.0 + dny = 0.0 + for i = 1:n + sk = atol[i] + rtol[i]*abs(y[i]) + dnf += (f0[i]/sk)^2 + dny += (y[i]/sk)^2 + end + if dnf <= 1e-10 || dny <= 1e-10 + h = 1e-6 + else + h = sqrt(dny/dnf)*0.01 + end + h = min(h, hmax) + h = h*posneg + y1 = y + h*y0 + f(f1, x+h, y1, params) + der2 = 0.0 + for i = 1:n + sk = atol[i] + rtol[i]*abs(y[i]) + der2 += ((f1[i]-f0[i])/sk)^2 + end + der2 = sqrt(der2)/h + der12 = max(abs(der2), sqrt(dnf)) + if der12 <= 1e-15 + h1 = max(1e-6, abs(h)*1e-3) + else + h1 = (0.01/der12)^(1.0/iord) + end + h = min(100*abs(h), h1, hmax) + return h*posneg +end + +function dopcore(m::Int64, f::Function, x::Float64, y::Vector{Float64}, h::Float64, k1::Vector{Float64}, k2::Vector{Float64}, k3::Vector{Float64}, k4::Vector{Float64}, k5::Vector{Float64}, k6::Vector{Float64}, k7::Vector{Float64}, k8::Vector{Float64}, k9::Vector{Float64}, k10::Vector{Float64}, atol::Vector{Float64}, rtol::Vector{Float64}, params::Any) + a21 = 5.26001519587677318785587544488e-2 + a31 = 1.97250569845378994544595329183e-2 + a32 = 5.91751709536136983633785987549e-2 + a41 = 2.95875854768068491816892993775e-2 + a43 = 8.87627564304205475450678981324e-2 + a51 = 2.41365134159266685502369798665e-1 + a53 = -8.84549479328286085344864962717e-1 + a54 = 9.24834003261792003115737966543e-1 + a61 = 3.7037037037037037037037037037e-2 + a64 = 1.70828608729473871279604482173e-1 + a65 = 1.25467687566822425016691814123e-1 + a71 = 3.7109375e-2 + a74 = 1.70252211019544039314978060272e-1 + a75 = 6.02165389804559606850219397283e-2 + a76 = -1.7578125e-2 + a81 = 3.70920001185047927108779319836e-2 + a84 = 1.70383925712239993810214054705e-1 + a85 = 1.07262030446373284651809199168e-1 + a86 = -1.53194377486244017527936158236e-2 + a87 = 8.27378916381402288758473766002e-3 + a91 = 6.24110958716075717114429577812e-1 + a94 = -3.36089262944694129406857109825e0 + a95 = -8.68219346841726006818189891453e-1 + a96 = 2.75920996994467083049415600797e1 + a97 = 2.01540675504778934086186788979e1 + a98 = -4.34898841810699588477366255144e1 + a101 = 4.77662536438264365890433908527e-1 + a104 = -2.48811461997166764192642586468e0 + a105 = -5.90290826836842996371446475743e-1 + a106 = 2.12300514481811942347288949897e1 + a107 = 1.52792336328824235832596922938e1 + a108 = -3.32882109689848629194453265587e1 + a109 = -2.03312017085086261358222928593e-2 + a111 = -9.3714243008598732571704021658e-1 + a114 = 5.18637242884406370830023853209e0 + a115 = 1.09143734899672957818500254654e0 + a116 = -8.14978701074692612513997267357e0 + a117 = -1.85200656599969598641566180701e1 + a118 = 2.27394870993505042818970056734e1 + a119 = 2.49360555267965238987089396762e0 + a1110 = -3.0467644718982195003823669022e0 + a121 = 2.27331014751653820792359768449e0 + a124 = -1.05344954667372501984066689879e1 + a125 = -2.00087205822486249909675718444e0 + a126 = -1.79589318631187989172765950534e1 + a127 = 2.79488845294199600508499808837e1 + a128 = -2.85899827713502369474065508674e0 + a129 = -8.87285693353062954433549289258e0 + a1210 = 1.23605671757943030647266201528e1 + a1211 = 6.43392746015763530355970484046e-1 + a141 = 5.61675022830479523392909219681e-2 + a147 = 2.53500210216624811088794765333e-1 + a148 = -2.46239037470802489917441475441e-1 + a149 = -1.24191423263816360469010140626e-1 + a1410 = 1.5329179827876569731206322685e-1 + a1411 = 8.20105229563468988491666602057e-3 + a1412 = 7.56789766054569976138603589584e-3 + a1413 = -8.298e-3 + a151 = 3.18346481635021405060768473261e-2 + a156 = 2.83009096723667755288322961402e-2 + a157 = 5.35419883074385676223797384372e-2 + a158 = -5.49237485713909884646569340306e-2 + a1511 = -1.08347328697249322858509316994e-4 + a1512 = 3.82571090835658412954920192323e-4 + a1513 = -3.40465008687404560802977114492e-4 + a1514 = 1.41312443674632500278074618366e-1 + a161 = -4.28896301583791923408573538692e-1 + a166 = -4.69762141536116384314449447206e0 + a167 = 7.68342119606259904184240953878e0 + a168 = 4.06898981839711007970213554331e0 + a169 = 3.56727187455281109270669543021e-1 + a1613 = -1.39902416515901462129418009734e-3 + a1614 = 2.9475147891527723389556272149e0 + a1615 = -9.15095847217987001081870187138e0 + b1 = 5.42937341165687622380535766363e-2 + b6 = 4.45031289275240888144113950566e0 + b7 = 1.89151789931450038304281599044e0 + b8 = -5.8012039600105847814672114227e0 + b9 = 3.1116436695781989440891606237e-1 + b10 = -1.52160949662516078556178806805e-1 + b11 = 2.01365400804030348374776537501e-1 + b12 = 4.47106157277725905176885569043e-2 + c2 = 0.526001519587677318785587544488e-01 + c3 = 0.789002279381515978178381316732e-01 + c4 = 0.118350341907227396726757197510e+00 + c5 = 0.281649658092772603273242802490e+00 + c6 = 0.333333333333333333333333333333e+00 + c7 = 0.25e+00 + c8 = 0.307692307692307692307692307692e+00 + c9 = 0.651282051282051282051282051282e+00 + c10 = 0.6e+00 + c11 = 0.857142857142857142857142857142e+00 + bhh1 = 0.244094488188976377952755905512e+00 + bhh2 = 0.733846688281611857341361741547e+00 + bhh3 = 0.220588235294117647058823529412e-01 + er1 = 0.1312004499419488073250102996e-01 + er6 = -0.1225156446376204440720569753e+01 + er7 = -0.4957589496572501915214079952e+00 + er8 = 0.1664377182454986536961530415e+01 + er9 = -0.3503288487499736816886487290e+00 + er10 = 0.3341791187130174790297318841e+00 + er11 = 0.8192320648511571246570742613e-01 + er12 = -0.2235530786388629525884427845e-01 + + #if m == 2 + # println(y) + #end + n=6 + y1 = zeros(n) + for i = 1:n + y1[i] = y[i] + h*a21*k1[i] + end + f(k2, x+c2*h, y1, params) + for i = 1:n + y1[i] = y[i] + h*(a31*k1[i] + a32*k2[i]) + end + f(k3, x+c3*h, y1, params) + for i = 1:n + y1[i] = y[i]+h*(a41*k1[i]+a43*k3[i]) + end + f(k4, x+c4*h, y1, params) + for i = 1:n + y1[i] = y[i]+h*(a51*k1[i]+a53*k3[i]+a54*k4[i]) + end + f(k5, x+c5*h, y1, params) + for i = 1:n + y1[i] = y[i]+h*(a61*k1[i]+a64*k4[i]+a65*k5[i]) + end + f(k6, x+c6*h, y1, params) + for i = 1:n + y1[i] = y[i]+h*(a71*k1[i]+a74*k4[i]+a75*k5[i]+a76*k6[i]) + end + f(k7, x+c7*h, y1, params) + for i = 1:n + y1[i] = y[i]+h*(a81*k1[i]+a84*k4[i]+a85*k5[i]+a86*k6[i]+a87*k7[i]) + end + f(k8, x+c8*h, y1, params) + for i = 1:n + y1[i] = y[i]+h*(a91*k1[i]+a94*k4[i]+a95*k5[i]+a96*k6[i]+a97*k7[i]+a98*k8[i]) + end + f(k9, x+c9*h, y1, params) + for i = 1:n + y1[i] = y[i]+h*(a101*k1[i]+a104*k4[i]+a105*k5[i]+a106*k6[i] + +a107*k7[i]+a108*k8[i]+a109*k9[i]) + end + f(k10, x+c10*h, y1, params) + for i = 1:n + y1[i] = y[i]+h*(a111*k1[i]+a114*k4[i]+a115*k5[i]+a116*k6[i] + +a117*k7[i]+a118*k8[i]+a119*k9[i]+a1110*k10[i]) + end + f(k2, x+c11*h, y1, params) + for i = 1:n + y1[i] = y[i]+h*(a121*k1[i]+a124*k4[i]+a125*k5[i]+a126*k6[i] + +a127*k7[i]+a128*k8[i]+a129*k9[i]+a1210*k10[i]+a1211*k2[i]) + end + f(k3, x+h, y1, params) + for i = 1:n + k4[i] = b1*k1[i]+b6*k6[i]+b7*k7[i]+b8*k8[i]+b9*k9[i]+b10*k10[i]+b11*k2[i]+b12*k3[i] + k5[i] = y[i]+h*k4[i] + end + + # Error estimation + err = 0.0 + err2 = 0.0 + for i = 1:n + sk = atol[i] + rtol[i]*max(abs(y[i]),abs(k5[i])) + erri = k4[i] - bhh1*k1[i] - bhh2*k9[i] - bhh3*k3[i] + err2 += (erri/sk)*(erri/sk) + erri = er1*k1[i] + er6*k6[i] + er7*k7[i] + er8*k8[i] + er9*k9[i] + er10*k10[i] + er11*k2[i] + er12*k3[i] + err += (erri/sk)*(erri/sk) + end + deno = err + 0.01*err2 + if deno <= 0.0 + deno = 1.0 + end + err = abs(h)*err*sqrt(1.0/(n*deno)) + return y1, err +end From 8a536f0f6b975bec0859eb2ea9fd29f5d5984c53 Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Mon, 24 Mar 2014 16:11:54 +0100 Subject: [PATCH 02/20] Remove debugging debris. --- src/dop853.jl | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/src/dop853.jl b/src/dop853.jl index 9309dcdf4..22faf9084 100644 --- a/src/dop853.jl +++ b/src/dop853.jl @@ -107,10 +107,6 @@ function dop853(f::Function, x::Float64, yin::Vector{Float64}, xend::Float64, pa end while true - #if nstep == 4 - # println(y) - # exit() - #end if nstep > nmax error("Exit at x=$x. More than nmax=$nmax steps needed.") end @@ -126,7 +122,7 @@ function dop853(f::Function, x::Float64, yin::Vector{Float64}, xend::Float64, pa f(k1, x, y, params) end - y1, err = dopcore(nstep, f, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, atol, rtol, params) + y1, err = dopcore(n, f, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, atol, rtol, params) xph = x+h nfcn += 11 @@ -237,7 +233,7 @@ function hinit(n::Int64, f::Function, x::Float64, y::Vector{Float64}, xend::Floa return h*posneg end -function dopcore(m::Int64, f::Function, x::Float64, y::Vector{Float64}, h::Float64, k1::Vector{Float64}, k2::Vector{Float64}, k3::Vector{Float64}, k4::Vector{Float64}, k5::Vector{Float64}, k6::Vector{Float64}, k7::Vector{Float64}, k8::Vector{Float64}, k9::Vector{Float64}, k10::Vector{Float64}, atol::Vector{Float64}, rtol::Vector{Float64}, params::Any) +function dopcore(n::Int64, f::Function, x::Float64, y::Vector{Float64}, h::Float64, k1::Vector{Float64}, k2::Vector{Float64}, k3::Vector{Float64}, k4::Vector{Float64}, k5::Vector{Float64}, k6::Vector{Float64}, k7::Vector{Float64}, k8::Vector{Float64}, k9::Vector{Float64}, k10::Vector{Float64}, atol::Vector{Float64}, rtol::Vector{Float64}, params::Any) a21 = 5.26001519587677318785587544488e-2 a31 = 1.97250569845378994544595329183e-2 a32 = 5.91751709536136983633785987549e-2 @@ -342,10 +338,6 @@ function dopcore(m::Int64, f::Function, x::Float64, y::Vector{Float64}, h::Float er11 = 0.8192320648511571246570742613e-01 er12 = -0.2235530786388629525884427845e-01 - #if m == 2 - # println(y) - #end - n=6 y1 = zeros(n) for i = 1:n y1[i] = y[i] + h*a21*k1[i] From 1bb01611fb754c90e8887b47eae0707a4f6b0ce8 Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Tue, 1 Apr 2014 11:06:24 +0200 Subject: [PATCH 03/20] Start API adjustment. --- src/dop853.jl | 71 +++++++++++++++++++++++++++++---------------------- 1 file changed, 40 insertions(+), 31 deletions(-) diff --git a/src/dop853.jl b/src/dop853.jl index 22faf9084..5e72afa23 100644 --- a/src/dop853.jl +++ b/src/dop853.jl @@ -1,8 +1,15 @@ -function dop853(f::Function, x::Float64, yin::Vector{Float64}, xend::Float64, params::Any; rtol::Vector{Float64}=[1e-6], atol::Vector{Float64}=[sqrt(eps())], - uround::Float64=2.3e-16, safe::Float64=0.9, fac1::Float64=0.333, fac2::Float64=6.0, - beta::Float64=0.0, maxstepsize::Float64=xend-x, initialstep::Float64=0.0, +abstract ODEProblem + +type ODEProblemFunction <: ODEProblem + f::Function +end + +function dop853(F::Function, y0, tspan; + reltol::Vector{Float64}=[1e-6], abstol::Vector{Float64}=[sqrt(eps())], + uround::Float64=eps(), safe::Float64=0.9, fac1::Float64=0.333, fac2::Float64=6.0, + beta::Float64=0.0, maxstepsize::Float64=tspan[end]-tspan[1], initialstep::Float64=0.0, maxsteps::Int64=100000, printmessages::Bool=false, nstiff::Int64=1000, - iout::Int64=0, solout::Function=s(x...)=return, dense::Vector{Int64}=[1:length(yin)]) + iout::Int64=0, solout::Function=s(x...)=return, dense::Vector{Int64}=[1:length(y0)]) c14 = 0.1e+00 c15 = 0.2e+00 @@ -58,7 +65,9 @@ function dop853(f::Function, x::Float64, yin::Vector{Float64}, xend::Float64, pa irtrn = 0 - n = length(yin) + x = tspan[1] + xend = tspan[end] + n = length(y0) nrd = length(dense) nfcn = 0 nstep = 0 @@ -75,32 +84,32 @@ function dop853(f::Function, x::Float64, yin::Vector{Float64}, xend::Float64, pa k9 = zeros(n) k10 = zeros(n) y = zeros(n) - copy!(y, yin) + copy!(y, y0) facold = 1e-4 expo1 = 1.0/8.0 - beta*0.2 facc1 = 1.0/fac1 facc2 = 1.0/fac2 posneg = sign(xend-x) - atol = length(atol) == 1 ? fill(atol[1], n) : atol - rtol = length(rtol) == 1 ? fill(rtol[1], n) : rtol + abstol = length(abstol) == 1 ? fill(abstol[1], n) : abstol + reltol = length(reltol) == 1 ? fill(reltol[1], n) : reltol last = false hlamb = 0.0 iasti = 0 - f(k1, x, y, params) + F(k1, x, y) nfcn += 1 hmax = abs(maxstepsize) nmax = maxsteps h = initialstep iord = 8 if h == 0.0 - h = hinit(n, f, x, y, xend, posneg, k1, k2, k3, iord, hmax, atol, rtol, params) + h = hinit(n, F, x, y, xend, posneg, k1, k2, k3, iord, hmax, abstol, reltol) nfcn += 1 end reject = false xold = x if iout != 0 #hout = 1.0 - irtrn = solout(naccpt+1, xold, x, y, con, icomp, params) + irtrn = solout(naccpt+1, xold, x, y, con, icomp) if irtrn < 0 return y end @@ -119,10 +128,10 @@ function dop853(f::Function, x::Float64, yin::Vector{Float64}, xend::Float64, pa end nstep += 1 if irtrn >= 2 - f(k1, x, y, params) + F(k1, x, y) end - y1, err = dopcore(n, f, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, atol, rtol, params) + y1, err = dopcore(n, F, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, abstol, reltol) xph = x+h nfcn += 11 @@ -133,7 +142,7 @@ function dop853(f::Function, x::Float64, yin::Vector{Float64}, xend::Float64, pa if err <= 1.0 facold = max(err, 1e-4) naccpt += 1 - f(k4, xph, k5, params) + F(k4, xph, k5,) nfcn += 1 # Stiffness detection if mod(naccpt, nstiff) == 0 || iasti > 0 @@ -200,11 +209,11 @@ function dop853(f::Function, x::Float64, yin::Vector{Float64}, xend::Float64, pa end end -function hinit(n::Int64, f::Function, x::Float64, y::Vector{Float64}, xend::Float64, posneg::Float64, f0::AbstractArray{Float64,1}, f1::AbstractArray{Float64,1}, y0::AbstractArray{Float64,1}, iord::Int64, hmax::Float64, atol::Vector{Float64}, rtol::Vector{Float64}, params::Any) +function hinit(n::Int64, F::Function, x::Float64, y::Vector{Float64}, xend::Float64, posneg::Float64, f0::AbstractArray{Float64,1}, f1::AbstractArray{Float64,1}, y0::AbstractArray{Float64,1}, iord::Int64, hmax::Float64, abstol::Vector{Float64}, reltol::Vector{Float64}) dnf = 0.0 dny = 0.0 for i = 1:n - sk = atol[i] + rtol[i]*abs(y[i]) + sk = abstol[i] + reltol[i]*abs(y[i]) dnf += (f0[i]/sk)^2 dny += (y[i]/sk)^2 end @@ -216,10 +225,10 @@ function hinit(n::Int64, f::Function, x::Float64, y::Vector{Float64}, xend::Floa h = min(h, hmax) h = h*posneg y1 = y + h*y0 - f(f1, x+h, y1, params) + F(f1, x+h, y1) der2 = 0.0 for i = 1:n - sk = atol[i] + rtol[i]*abs(y[i]) + sk = abstol[i] + reltol[i]*abs(y[i]) der2 += ((f1[i]-f0[i])/sk)^2 end der2 = sqrt(der2)/h @@ -233,7 +242,7 @@ function hinit(n::Int64, f::Function, x::Float64, y::Vector{Float64}, xend::Floa return h*posneg end -function dopcore(n::Int64, f::Function, x::Float64, y::Vector{Float64}, h::Float64, k1::Vector{Float64}, k2::Vector{Float64}, k3::Vector{Float64}, k4::Vector{Float64}, k5::Vector{Float64}, k6::Vector{Float64}, k7::Vector{Float64}, k8::Vector{Float64}, k9::Vector{Float64}, k10::Vector{Float64}, atol::Vector{Float64}, rtol::Vector{Float64}, params::Any) +function dopcore(n::Int64, F::Function, x::Float64, y::Vector{Float64}, h::Float64, k1::Vector{Float64}, k2::Vector{Float64}, k3::Vector{Float64}, k4::Vector{Float64}, k5::Vector{Float64}, k6::Vector{Float64}, k7::Vector{Float64}, k8::Vector{Float64}, k9::Vector{Float64}, k10::Vector{Float64}, abstol::Vector{Float64}, reltol::Vector{Float64}) a21 = 5.26001519587677318785587544488e-2 a31 = 1.97250569845378994544595329183e-2 a32 = 5.91751709536136983633785987549e-2 @@ -342,50 +351,50 @@ function dopcore(n::Int64, f::Function, x::Float64, y::Vector{Float64}, h::Float for i = 1:n y1[i] = y[i] + h*a21*k1[i] end - f(k2, x+c2*h, y1, params) + F(k2, x+c2*h, y1) for i = 1:n y1[i] = y[i] + h*(a31*k1[i] + a32*k2[i]) end - f(k3, x+c3*h, y1, params) + F(k3, x+c3*h, y1) for i = 1:n y1[i] = y[i]+h*(a41*k1[i]+a43*k3[i]) end - f(k4, x+c4*h, y1, params) + F(k4, x+c4*h, y1) for i = 1:n y1[i] = y[i]+h*(a51*k1[i]+a53*k3[i]+a54*k4[i]) end - f(k5, x+c5*h, y1, params) + F(k5, x+c5*h, y1) for i = 1:n y1[i] = y[i]+h*(a61*k1[i]+a64*k4[i]+a65*k5[i]) end - f(k6, x+c6*h, y1, params) + F(k6, x+c6*h, y1) for i = 1:n y1[i] = y[i]+h*(a71*k1[i]+a74*k4[i]+a75*k5[i]+a76*k6[i]) end - f(k7, x+c7*h, y1, params) + F(k7, x+c7*h, y1) for i = 1:n y1[i] = y[i]+h*(a81*k1[i]+a84*k4[i]+a85*k5[i]+a86*k6[i]+a87*k7[i]) end - f(k8, x+c8*h, y1, params) + F(k8, x+c8*h, y1) for i = 1:n y1[i] = y[i]+h*(a91*k1[i]+a94*k4[i]+a95*k5[i]+a96*k6[i]+a97*k7[i]+a98*k8[i]) end - f(k9, x+c9*h, y1, params) + F(k9, x+c9*h, y1) for i = 1:n y1[i] = y[i]+h*(a101*k1[i]+a104*k4[i]+a105*k5[i]+a106*k6[i] +a107*k7[i]+a108*k8[i]+a109*k9[i]) end - f(k10, x+c10*h, y1, params) + F(k10, x+c10*h, y1) for i = 1:n y1[i] = y[i]+h*(a111*k1[i]+a114*k4[i]+a115*k5[i]+a116*k6[i] +a117*k7[i]+a118*k8[i]+a119*k9[i]+a1110*k10[i]) end - f(k2, x+c11*h, y1, params) + F(k2, x+c11*h, y1) for i = 1:n y1[i] = y[i]+h*(a121*k1[i]+a124*k4[i]+a125*k5[i]+a126*k6[i] +a127*k7[i]+a128*k8[i]+a129*k9[i]+a1210*k10[i]+a1211*k2[i]) end - f(k3, x+h, y1, params) + F(k3, x+h, y1) for i = 1:n k4[i] = b1*k1[i]+b6*k6[i]+b7*k7[i]+b8*k8[i]+b9*k9[i]+b10*k10[i]+b11*k2[i]+b12*k3[i] k5[i] = y[i]+h*k4[i] @@ -395,7 +404,7 @@ function dopcore(n::Int64, f::Function, x::Float64, y::Vector{Float64}, h::Float err = 0.0 err2 = 0.0 for i = 1:n - sk = atol[i] + rtol[i]*max(abs(y[i]),abs(k5[i])) + sk = abstol[i] + reltol[i]*max(abs(y[i]),abs(k5[i])) erri = k4[i] - bhh1*k1[i] - bhh2*k9[i] - bhh3*k3[i] err2 += (erri/sk)*(erri/sk) erri = er1*k1[i] + er6*k6[i] + er7*k7[i] + er8*k8[i] + er9*k9[i] + er10*k10[i] + er11*k2[i] + er12*k3[i] From 331f2b08dc10be41db3c7eda80717053d9739082 Mon Sep 17 00:00:00 2001 From: Andrei Berceanu Date: Mon, 23 Jun 2014 15:43:34 +0200 Subject: [PATCH 04/20] initial commit --- src/dop853.jl | 66 +++++++++++++++++++++++++++------------------------ 1 file changed, 35 insertions(+), 31 deletions(-) diff --git a/src/dop853.jl b/src/dop853.jl index 5e72afa23..54c95d406 100644 --- a/src/dop853.jl +++ b/src/dop853.jl @@ -7,7 +7,7 @@ end function dop853(F::Function, y0, tspan; reltol::Vector{Float64}=[1e-6], abstol::Vector{Float64}=[sqrt(eps())], uround::Float64=eps(), safe::Float64=0.9, fac1::Float64=0.333, fac2::Float64=6.0, - beta::Float64=0.0, maxstepsize::Float64=tspan[end]-tspan[1], initialstep::Float64=0.0, + beta::Float64=0.0, maxstepsize::Float64=tspan[end]-tspan[1], initialstep::Float64=1e-7, maxsteps::Int64=100000, printmessages::Bool=false, nstiff::Int64=1000, iout::Int64=0, solout::Function=s(x...)=return, dense::Vector{Int64}=[1:length(y0)]) @@ -67,23 +67,23 @@ function dop853(F::Function, y0, tspan; x = tspan[1] xend = tspan[end] - n = length(y0) + n = length(y0) #1 nrd = length(dense) nfcn = 0 nstep = 0 naccpt = 0 nrejct = 0 - k1 = zeros(n) - k2 = zeros(n) - k3 = zeros(n) - k4 = zeros(n) - k5 = zeros(n) - k6 = zeros(n) - k7 = zeros(n) - k8 = zeros(n) - k9 = zeros(n) - k10 = zeros(n) - y = zeros(n) + k1 = 0. * y0 + k2 = 0. * y0 + k3 = 0. * y0 + k4 = 0. * y0 + k5 = 0. * y0 + k6 = 0. * y0 + k7 = 0. * y0 + k8 = 0. * y0 + k9 = 0. * y0 + k10 = 0. * y0 + y = 0. * y0 copy!(y, y0) facold = 1e-4 expo1 = 1.0/8.0 - beta*0.2 @@ -103,6 +103,7 @@ function dop853(F::Function, y0, tspan; iord = 8 if h == 0.0 h = hinit(n, F, x, y, xend, posneg, k1, k2, k3, iord, hmax, abstol, reltol) + println("hinit = $h") nfcn += 1 end reject = false @@ -132,6 +133,7 @@ function dop853(F::Function, y0, tspan; end y1, err = dopcore(n, F, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, abstol, reltol) + println("error is $err") xph = x+h nfcn += 11 @@ -142,7 +144,7 @@ function dop853(F::Function, y0, tspan; if err <= 1.0 facold = max(err, 1e-4) naccpt += 1 - F(k4, xph, k5,) + F(k4, xph, k5) nfcn += 1 # Stiffness detection if mod(naccpt, nstiff) == 0 || iasti > 0 @@ -181,6 +183,7 @@ function dop853(F::Function, y0, tspan; copy!(y, k5) xold = x x = xph + println(x) # if # solout # end @@ -191,12 +194,12 @@ function dop853(F::Function, y0, tspan; return y end if abs(hnew) > hmax - hnew = posneg*hmax + hnew = posneg*hmax end - if reject + if reject hnew = posneg*min(abs(hnew),abs(h)) end - reject = false + reject = false else hnew = h/min(facc1,fac11/safe) reject = true @@ -209,7 +212,8 @@ function dop853(F::Function, y0, tspan; end end -function hinit(n::Int64, F::Function, x::Float64, y::Vector{Float64}, xend::Float64, posneg::Float64, f0::AbstractArray{Float64,1}, f1::AbstractArray{Float64,1}, y0::AbstractArray{Float64,1}, iord::Int64, hmax::Float64, abstol::Vector{Float64}, reltol::Vector{Float64}) +#f0, f1, y0 -> k1, k2, k3 +function hinit(n::Int64, F::Function, x::Float64, y::Vector, xend::Float64, posneg::Float64, f0::AbstractVector, f1::AbstractVector, y0::AbstractVector, iord::Int64, hmax::Float64, abstol::Vector{Float64}, reltol::Vector{Float64}) dnf = 0.0 dny = 0.0 for i = 1:n @@ -217,10 +221,10 @@ function hinit(n::Int64, F::Function, x::Float64, y::Vector{Float64}, xend::Floa dnf += (f0[i]/sk)^2 dny += (y[i]/sk)^2 end - if dnf <= 1e-10 || dny <= 1e-10 + if maximum(abs(dnf)) <= 1e-10 || maximum(abs(dny)) <= 1e-10 h = 1e-6 else - h = sqrt(dny/dnf)*0.01 + h = sqrt(maximum(abs(dny/dnf)))*0.01 end h = min(h, hmax) h = h*posneg @@ -231,8 +235,8 @@ function hinit(n::Int64, F::Function, x::Float64, y::Vector{Float64}, xend::Floa sk = abstol[i] + reltol[i]*abs(y[i]) der2 += ((f1[i]-f0[i])/sk)^2 end - der2 = sqrt(der2)/h - der12 = max(abs(der2), sqrt(dnf)) + der2 = sqrt(maximum(abs(der2)))/h + der12 = max(maximum(abs(der2)), sqrt(maximum(abs(dnf)))) if der12 <= 1e-15 h1 = max(1e-6, abs(h)*1e-3) else @@ -242,7 +246,7 @@ function hinit(n::Int64, F::Function, x::Float64, y::Vector{Float64}, xend::Floa return h*posneg end -function dopcore(n::Int64, F::Function, x::Float64, y::Vector{Float64}, h::Float64, k1::Vector{Float64}, k2::Vector{Float64}, k3::Vector{Float64}, k4::Vector{Float64}, k5::Vector{Float64}, k6::Vector{Float64}, k7::Vector{Float64}, k8::Vector{Float64}, k9::Vector{Float64}, k10::Vector{Float64}, abstol::Vector{Float64}, reltol::Vector{Float64}) +function dopcore(n::Int64, F::Function, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector{Float64}, reltol::Vector{Float64}) a21 = 5.26001519587677318785587544488e-2 a31 = 1.97250569845378994544595329183e-2 a32 = 5.91751709536136983633785987549e-2 @@ -346,8 +350,8 @@ function dopcore(n::Int64, F::Function, x::Float64, y::Vector{Float64}, h::Float er10 = 0.3341791187130174790297318841e+00 er11 = 0.8192320648511571246570742613e-01 er12 = -0.2235530786388629525884427845e-01 - - y1 = zeros(n) + + y1 = 0. * y for i = 1:n y1[i] = y[i] + h*a21*k1[i] end @@ -357,7 +361,7 @@ function dopcore(n::Int64, F::Function, x::Float64, y::Vector{Float64}, h::Float end F(k3, x+c3*h, y1) for i = 1:n - y1[i] = y[i]+h*(a41*k1[i]+a43*k3[i]) + y1[i] = y[i]+h*(a41*k1[i]+a43*k3[i]) end F(k4, x+c4*h, y1) for i = 1:n @@ -373,7 +377,7 @@ function dopcore(n::Int64, F::Function, x::Float64, y::Vector{Float64}, h::Float end F(k7, x+c7*h, y1) for i = 1:n - y1[i] = y[i]+h*(a81*k1[i]+a84*k4[i]+a85*k5[i]+a86*k6[i]+a87*k7[i]) + y1[i] = y[i]+h*(a81*k1[i]+a84*k4[i]+a85*k5[i]+a86*k6[i]+a87*k7[i]) end F(k8, x+c8*h, y1) for i = 1:n @@ -404,16 +408,16 @@ function dopcore(n::Int64, F::Function, x::Float64, y::Vector{Float64}, h::Float err = 0.0 err2 = 0.0 for i = 1:n - sk = abstol[i] + reltol[i]*max(abs(y[i]),abs(k5[i])) + sk = abstol[i] + reltol[i]*max(maximum(abs(y[i])), maximum(abs(k5[i]))) erri = k4[i] - bhh1*k1[i] - bhh2*k9[i] - bhh3*k3[i] err2 += (erri/sk)*(erri/sk) erri = er1*k1[i] + er6*k6[i] + er7*k7[i] + er8*k8[i] + er9*k9[i] + er10*k10[i] + er11*k2[i] + er12*k3[i] err += (erri/sk)*(erri/sk) - end - deno = err + 0.01*err2 + end + deno = maximum(abs(err)) + 0.01*maximum(abs(err2)) if deno <= 0.0 deno = 1.0 end err = abs(h)*err*sqrt(1.0/(n*deno)) - return y1, err + return y1, maximum(abs(err)) end From b9ba1849710a189c953c5b380d9ecbfa1a17b2f2 Mon Sep 17 00:00:00 2001 From: Andrei Berceanu Date: Mon, 23 Jun 2014 15:56:01 +0200 Subject: [PATCH 05/20] reset initial step to default value --- src/dop853.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dop853.jl b/src/dop853.jl index 54c95d406..30ec44033 100644 --- a/src/dop853.jl +++ b/src/dop853.jl @@ -7,7 +7,7 @@ end function dop853(F::Function, y0, tspan; reltol::Vector{Float64}=[1e-6], abstol::Vector{Float64}=[sqrt(eps())], uround::Float64=eps(), safe::Float64=0.9, fac1::Float64=0.333, fac2::Float64=6.0, - beta::Float64=0.0, maxstepsize::Float64=tspan[end]-tspan[1], initialstep::Float64=1e-7, + beta::Float64=0.0, maxstepsize::Float64=tspan[end]-tspan[1], initialstep::Float64=0.0, maxsteps::Int64=100000, printmessages::Bool=false, nstiff::Int64=1000, iout::Int64=0, solout::Function=s(x...)=return, dense::Vector{Int64}=[1:length(y0)]) From f2c6fa4ced6351a0bbd0dfbad1d9b37f14a7cd20 Mon Sep 17 00:00:00 2001 From: Andrei Berceanu Date: Mon, 23 Jun 2014 16:15:08 +0200 Subject: [PATCH 06/20] printing debug messages is now optional --- src/dop853.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/dop853.jl b/src/dop853.jl index 30ec44033..137dfcfec 100644 --- a/src/dop853.jl +++ b/src/dop853.jl @@ -103,7 +103,7 @@ function dop853(F::Function, y0, tspan; iord = 8 if h == 0.0 h = hinit(n, F, x, y, xend, posneg, k1, k2, k3, iord, hmax, abstol, reltol) - println("hinit = $h") + printmessages && println("hinit = $h") nfcn += 1 end reject = false @@ -133,7 +133,7 @@ function dop853(F::Function, y0, tspan; end y1, err = dopcore(n, F, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, abstol, reltol) - println("error is $err") + printmessages && println("error is $err") xph = x+h nfcn += 11 @@ -183,7 +183,7 @@ function dop853(F::Function, y0, tspan; copy!(y, k5) xold = x x = xph - println(x) + printmessages && println(x) # if # solout # end From bc07775a45627b6d173ff048b810874f4ae837e2 Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Tue, 24 Jun 2014 15:01:32 +0200 Subject: [PATCH 07/20] API adjustments, generalizations and performance tweaking. --- src/dop853.jl | 49 +++++++++++++++++++++++++++++++++++-------------- 1 file changed, 35 insertions(+), 14 deletions(-) diff --git a/src/dop853.jl b/src/dop853.jl index 5e72afa23..1a3fd7f46 100644 --- a/src/dop853.jl +++ b/src/dop853.jl @@ -4,12 +4,13 @@ type ODEProblemFunction <: ODEProblem f::Function end -function dop853(F::Function, y0, tspan; - reltol::Vector{Float64}=[1e-6], abstol::Vector{Float64}=[sqrt(eps())], - uround::Float64=eps(), safe::Float64=0.9, fac1::Float64=0.333, fac2::Float64=6.0, - beta::Float64=0.0, maxstepsize::Float64=tspan[end]-tspan[1], initialstep::Float64=0.0, - maxsteps::Int64=100000, printmessages::Bool=false, nstiff::Int64=1000, - iout::Int64=0, solout::Function=s(x...)=return, dense::Vector{Int64}=[1:length(y0)]) +function dop853(F, y0, tspan; + reltol=[1e-6], abstol=[sqrt(eps())], + safe=0.9, fac1=0.333, fac2=6.0, beta=0.0, + maxstep=tspan[end]-tspan[1], initstep=0.0, + maxsteps=100000, printmessages=false, nstiff=1000, + iout=0, solout=s(x...)=return, dense=[1:length(y0)], + points=:all) c14 = 0.1e+00 c15 = 0.2e+00 @@ -67,6 +68,9 @@ function dop853(F::Function, y0, tspan; x = tspan[1] xend = tspan[end] + tout = Array(typeof(tspan[1]),0) + yout = Array(typeof(y0*one(x)),0) + n = length(y0) nrd = length(dense) nfcn = 0 @@ -84,6 +88,7 @@ function dop853(F::Function, y0, tspan; k9 = zeros(n) k10 = zeros(n) y = zeros(n) + y1 = zeros(n) copy!(y, y0) facold = 1e-4 expo1 = 1.0/8.0 - beta*0.2 @@ -97,9 +102,9 @@ function dop853(F::Function, y0, tspan; iasti = 0 F(k1, x, y) nfcn += 1 - hmax = abs(maxstepsize) + hmax = abs(maxstep) nmax = maxsteps - h = initialstep + h = initstep iord = 8 if h == 0.0 h = hinit(n, F, x, y, xend, posneg, k1, k2, k3, iord, hmax, abstol, reltol) @@ -119,7 +124,7 @@ function dop853(F::Function, y0, tspan; if nstep > nmax error("Exit at x=$x. More than nmax=$nmax steps needed.") end - if 0.1*abs(h) <= abs(x)*uround + if 0.1*abs(h) <= abs(x)*eps() error("Exit at x=$x. Step size too small, h=$h.") end if (x+1.01*h-xend)*posneg > 0.0 @@ -131,7 +136,7 @@ function dop853(F::Function, y0, tspan; F(k1, x, y) end - y1, err = dopcore(n, F, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, abstol, reltol) + err = dopcore(y1, n, F, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, abstol, reltol) xph = x+h nfcn += 11 @@ -181,6 +186,10 @@ function dop853(F::Function, y0, tspan; copy!(y, k5) xold = x x = xph + if points == :all + push!(tout, x) + push!(yout, y1) + end # if # solout # end @@ -188,7 +197,11 @@ function dop853(F::Function, y0, tspan; if last h = hnew idid = 1 - return y + if points == :last + push!(yout, y1) + push!(tout, x) + end + return yout, tout end if abs(hnew) > hmax hnew = posneg*hmax @@ -209,6 +222,15 @@ function dop853(F::Function, y0, tspan; end end +function denseout(ind, t, told, h, coeff, dense) + if ~in(ind, dense) + error("No dense output available for component: $ind.") + else + s = (t - told)/h + s1 = 1.0 - s + end +end + function hinit(n::Int64, F::Function, x::Float64, y::Vector{Float64}, xend::Float64, posneg::Float64, f0::AbstractArray{Float64,1}, f1::AbstractArray{Float64,1}, y0::AbstractArray{Float64,1}, iord::Int64, hmax::Float64, abstol::Vector{Float64}, reltol::Vector{Float64}) dnf = 0.0 dny = 0.0 @@ -242,7 +264,7 @@ function hinit(n::Int64, F::Function, x::Float64, y::Vector{Float64}, xend::Floa return h*posneg end -function dopcore(n::Int64, F::Function, x::Float64, y::Vector{Float64}, h::Float64, k1::Vector{Float64}, k2::Vector{Float64}, k3::Vector{Float64}, k4::Vector{Float64}, k5::Vector{Float64}, k6::Vector{Float64}, k7::Vector{Float64}, k8::Vector{Float64}, k9::Vector{Float64}, k10::Vector{Float64}, abstol::Vector{Float64}, reltol::Vector{Float64}) +function dopcore(y1::Vector, n::Int64, F::Function, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector, reltol::Vector) a21 = 5.26001519587677318785587544488e-2 a31 = 1.97250569845378994544595329183e-2 a32 = 5.91751709536136983633785987549e-2 @@ -347,7 +369,6 @@ function dopcore(n::Int64, F::Function, x::Float64, y::Vector{Float64}, h::Float er11 = 0.8192320648511571246570742613e-01 er12 = -0.2235530786388629525884427845e-01 - y1 = zeros(n) for i = 1:n y1[i] = y[i] + h*a21*k1[i] end @@ -415,5 +436,5 @@ function dopcore(n::Int64, F::Function, x::Float64, y::Vector{Float64}, h::Float deno = 1.0 end err = abs(h)*err*sqrt(1.0/(n*deno)) - return y1, err + return err end From ae12a42143e938624186cb55b1aa75aa6290ab32 Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Wed, 25 Jun 2014 15:42:27 +0200 Subject: [PATCH 08/20] Implement new high-level API. --- src/ODE.jl | 2 +- src/dop853.jl | 46 +++++++++++++++++++++++++--------------------- 2 files changed, 26 insertions(+), 22 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 828386f88..59e50f341 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -5,7 +5,7 @@ module ODE using Polynomial ## minimal function export list -export ode23, ode4, ode45, ode4s, ode4ms, dop853 +export ode23, ode4, ode45, ode4s, ode4ms, dop853, ODEProblem ## complete function export list #export ode23, ode4, diff --git a/src/dop853.jl b/src/dop853.jl index 4ccf0de36..4a3eb87ca 100644 --- a/src/dop853.jl +++ b/src/dop853.jl @@ -4,7 +4,11 @@ type ODEProblemFunction <: ODEProblem f::Function end -function dop853(F, y0, tspan; +F!(p::ODEProblemFunction, y, x, t) = copy!(y, p.f(x, t)) + +dop853(f::Function, y0, tspan; args...) = dop853(F!, ODEProblemFunction(f), y0, tspan; args...) + +function dop853(F!::Function, p::ODEProblem, y0, tspan; reltol=[1e-6], abstol=[sqrt(eps())], safe=0.9, fac1=0.333, fac2=6.0, beta=0.0, maxstep=tspan[end]-tspan[1], initstep=0.0, @@ -101,14 +105,14 @@ function dop853(F, y0, tspan; last = false hlamb = 0.0 iasti = 0 - F(k1, x, y) + F!(p, k1, x, y) nfcn += 1 hmax = abs(maxstep) nmax = maxsteps h = initstep iord = 8 if h == 0.0 - h = hinit(n, F, x, y, xend, posneg, k1, k2, k3, iord, hmax, abstol, reltol) + h = hinit(n, F!, p, x, y, xend, posneg, k1, k2, k3, iord, hmax, abstol, reltol) printmessages && println("hinit = $h") nfcn += 1 end @@ -135,10 +139,10 @@ function dop853(F, y0, tspan; end nstep += 1 if irtrn >= 2 - F(k1, x, y) + F!(p, k1, x, y) end - err = dopcore(y1, n, F, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, abstol, reltol) + err = dopcore(y1, n, F!, p, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, abstol, reltol) xph = x+h nfcn += 11 @@ -149,7 +153,7 @@ function dop853(F, y0, tspan; if err <= 1.0 facold = max(err, 1e-4) naccpt += 1 - F(k4, xph, k5) + F!(p, k4, xph, k5) nfcn += 1 # Stiffness detection if mod(naccpt, nstiff) == 0 || iasti > 0 @@ -233,7 +237,7 @@ function denseout(ind, t, told, h, coeff, dense) end end -function hinit(n::Int64, F::Function, x::Float64, y::Vector, xend::Float64, posneg::Float64, f0::AbstractVector, f1::AbstractVector, y0::AbstractVector, iord::Int64, hmax::Float64, abstol::Vector{Float64}, reltol::Vector{Float64}) +function hinit(n::Int64, F!::Function, p::ODEProblem, x::Float64, y::Vector, xend::Float64, posneg::Float64, f0::AbstractVector, f1::AbstractVector, y0::AbstractVector, iord::Int64, hmax::Float64, abstol::Vector{Float64}, reltol::Vector{Float64}) dnf = 0.0 dny = 0.0 for i = 1:n @@ -249,7 +253,7 @@ function hinit(n::Int64, F::Function, x::Float64, y::Vector, xend::Float64, posn h = min(h, hmax) h = h*posneg y1 = y + h*y0 - F(f1, x+h, y1) + F!(p, f1, x+h, y1) der2 = 0.0 for i = 1:n sk = abstol[i] + reltol[i]*abs(y[i]) @@ -266,7 +270,7 @@ function hinit(n::Int64, F::Function, x::Float64, y::Vector, xend::Float64, posn return h*posneg end -function dopcore(y1::Vector, n::Int64, F::Function, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector, reltol::Vector) +function dopcore(y1::Vector, n::Int64, F!::Function, p::ODEProblem, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector, reltol::Vector) a21 = 5.26001519587677318785587544488e-2 a31 = 1.97250569845378994544595329183e-2 a32 = 5.91751709536136983633785987549e-2 @@ -370,54 +374,54 @@ function dopcore(y1::Vector, n::Int64, F::Function, x::Float64, y::Vector, h::Fl er10 = 0.3341791187130174790297318841e+00 er11 = 0.8192320648511571246570742613e-01 er12 = -0.2235530786388629525884427845e-01 - + for i = 1:n y1[i] = y[i] + h*a21*k1[i] end - F(k2, x+c2*h, y1) + F!(p, k2, x+c2*h, y1) for i = 1:n y1[i] = y[i] + h*(a31*k1[i] + a32*k2[i]) end - F(k3, x+c3*h, y1) + F!(p, k3, x+c3*h, y1) for i = 1:n y1[i] = y[i]+h*(a41*k1[i]+a43*k3[i]) end - F(k4, x+c4*h, y1) + F!(p, k4, x+c4*h, y1) for i = 1:n y1[i] = y[i]+h*(a51*k1[i]+a53*k3[i]+a54*k4[i]) end - F(k5, x+c5*h, y1) + F!(p, k5, x+c5*h, y1) for i = 1:n y1[i] = y[i]+h*(a61*k1[i]+a64*k4[i]+a65*k5[i]) end - F(k6, x+c6*h, y1) + F!(p, k6, x+c6*h, y1) for i = 1:n y1[i] = y[i]+h*(a71*k1[i]+a74*k4[i]+a75*k5[i]+a76*k6[i]) end - F(k7, x+c7*h, y1) + F!(p, k7, x+c7*h, y1) for i = 1:n y1[i] = y[i]+h*(a81*k1[i]+a84*k4[i]+a85*k5[i]+a86*k6[i]+a87*k7[i]) end - F(k8, x+c8*h, y1) + F!(p, k8, x+c8*h, y1) for i = 1:n y1[i] = y[i]+h*(a91*k1[i]+a94*k4[i]+a95*k5[i]+a96*k6[i]+a97*k7[i]+a98*k8[i]) end - F(k9, x+c9*h, y1) + F!(p, k9, x+c9*h, y1) for i = 1:n y1[i] = y[i]+h*(a101*k1[i]+a104*k4[i]+a105*k5[i]+a106*k6[i] +a107*k7[i]+a108*k8[i]+a109*k9[i]) end - F(k10, x+c10*h, y1) + F!(p, k10, x+c10*h, y1) for i = 1:n y1[i] = y[i]+h*(a111*k1[i]+a114*k4[i]+a115*k5[i]+a116*k6[i] +a117*k7[i]+a118*k8[i]+a119*k9[i]+a1110*k10[i]) end - F(k2, x+c11*h, y1) + F!(p, k2, x+c11*h, y1) for i = 1:n y1[i] = y[i]+h*(a121*k1[i]+a124*k4[i]+a125*k5[i]+a126*k6[i] +a127*k7[i]+a128*k8[i]+a129*k9[i]+a1210*k10[i]+a1211*k2[i]) end - F(k3, x+h, y1) + F!(p, k3, x+h, y1) for i = 1:n k4[i] = b1*k1[i]+b6*k6[i]+b7*k7[i]+b8*k8[i]+b9*k9[i]+b10*k10[i]+b11*k2[i]+b12*k3[i] k5[i] = y[i]+h*k4[i] From cd6ae6fb5f8f934b1df8115c513070e833ad301d Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Wed, 25 Jun 2014 17:09:05 +0200 Subject: [PATCH 09/20] Test suite integration and minor fixes. --- src/dop853.jl | 29 ++++++++++++++++++++--------- test/runtests.jl | 4 +++- 2 files changed, 23 insertions(+), 10 deletions(-) diff --git a/src/dop853.jl b/src/dop853.jl index 4a3eb87ca..99736fb10 100644 --- a/src/dop853.jl +++ b/src/dop853.jl @@ -4,9 +4,14 @@ type ODEProblemFunction <: ODEProblem f::Function end -F!(p::ODEProblemFunction, y, x, t) = copy!(y, p.f(x, t)) +F!(p::ODEProblemFunction, y, x, t) = copy!(y, [p.f(x, t)]) -dop853(f::Function, y0, tspan; args...) = dop853(F!, ODEProblemFunction(f), y0, tspan; args...) +dop853(f::Function, y0::Vector, tspan; args...) = dop853(F!, ODEProblemFunction(f), y0, tspan; args...) + +function dop853(f::Function, y0::Number, tspan; args...) + tout, yout = dop853(F!, ODEProblemFunction(f), [y0], tspan; args...) + return tout, vcat(yout...) +end function dop853(F!::Function, p::ODEProblem, y0, tspan; reltol=[1e-6], abstol=[sqrt(eps())], @@ -16,6 +21,7 @@ function dop853(F!::Function, p::ODEProblem, y0, tspan; iout=0, solout=s(x...)=return, dense=[1:length(y0)], points=:all) + c14 = 0.1e+00 c15 = 0.2e+00 c16 = 0.777777777777777777777777777778e+00 @@ -75,6 +81,11 @@ function dop853(F!::Function, p::ODEProblem, y0, tspan; tout = Array(typeof(tspan[1]),0) yout = Array(typeof(y0*one(x)),0) + if points == :all + push!(tout, x) + push!(yout, y0) + end + n = length(y0) nrd = length(dense) nfcn = 0 @@ -92,7 +103,6 @@ function dop853(F!::Function, p::ODEProblem, y0, tspan; k9 = 0. * y0 k10 = 0. * y0 y = 0. * y0 - y1 = 0. * y0 copy!(y, y0) facold = 1e-4 @@ -142,7 +152,7 @@ function dop853(F!::Function, p::ODEProblem, y0, tspan; F!(p, k1, x, y) end - err = dopcore(y1, n, F!, p, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, abstol, reltol) + y1, err = dopcore(n, F!, p, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, abstol, reltol) xph = x+h nfcn += 11 @@ -204,10 +214,9 @@ function dop853(F!::Function, p::ODEProblem, y0, tspan; h = hnew idid = 1 if points == :last - push!(yout, y1) - push!(tout, x) + return x, y1 end - return yout, tout + return tout, yout end if abs(hnew) > hmax hnew = posneg*hmax @@ -270,7 +279,7 @@ function hinit(n::Int64, F!::Function, p::ODEProblem, x::Float64, y::Vector, xen return h*posneg end -function dopcore(y1::Vector, n::Int64, F!::Function, p::ODEProblem, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector, reltol::Vector) +function dopcore(n::Int64, F!::Function, p::ODEProblem, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector, reltol::Vector) a21 = 5.26001519587677318785587544488e-2 a31 = 1.97250569845378994544595329183e-2 a32 = 5.91751709536136983633785987549e-2 @@ -375,6 +384,8 @@ function dopcore(y1::Vector, n::Int64, F!::Function, p::ODEProblem, x::Float64, er11 = 0.8192320648511571246570742613e-01 er12 = -0.2235530786388629525884427845e-01 + y1 = 0.0 * y + for i = 1:n y1[i] = y[i] + h*a21*k1[i] end @@ -442,5 +453,5 @@ function dopcore(y1::Vector, n::Int64, F!::Function, p::ODEProblem, x::Float64, deno = 1.0 end err = abs(h)*err*sqrt(1.0/(n*deno)) - return maximum(abs(err)) + return y1, maximum(abs(err)) end diff --git a/test/runtests.jl b/test/runtests.jl index def8d1774..2d55a401b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,7 +15,9 @@ solvers = [ ODE.ode4ms, ODE.ode4s_s, ODE.ode4s_kr, - ODE.ode78_fb] + ODE.ode78_fb, + + ODE.dop853] for solver in solvers println("using $solver") From 84cdf77c4959097c9b14bc7c7d0b5f66304c5868 Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Wed, 25 Jun 2014 17:37:21 +0200 Subject: [PATCH 10/20] Update low-level API. --- src/dop853.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/dop853.jl b/src/dop853.jl index 99736fb10..77d905707 100644 --- a/src/dop853.jl +++ b/src/dop853.jl @@ -6,14 +6,14 @@ end F!(p::ODEProblemFunction, y, x, t) = copy!(y, [p.f(x, t)]) -dop853(f::Function, y0::Vector, tspan; args...) = dop853(F!, ODEProblemFunction(f), y0, tspan; args...) +dop853(f::Function, y0::Vector, tspan; args...) = dop853(ODEProblemFunction(f), y0, tspan; args...) function dop853(f::Function, y0::Number, tspan; args...) tout, yout = dop853(F!, ODEProblemFunction(f), [y0], tspan; args...) return tout, vcat(yout...) end -function dop853(F!::Function, p::ODEProblem, y0, tspan; +function dop853(p::ODEProblem, y0, tspan; reltol=[1e-6], abstol=[sqrt(eps())], safe=0.9, fac1=0.333, fac2=6.0, beta=0.0, maxstep=tspan[end]-tspan[1], initstep=0.0, @@ -122,7 +122,7 @@ function dop853(F!::Function, p::ODEProblem, y0, tspan; h = initstep iord = 8 if h == 0.0 - h = hinit(n, F!, p, x, y, xend, posneg, k1, k2, k3, iord, hmax, abstol, reltol) + h = hinit(n, p, x, y, xend, posneg, k1, k2, k3, iord, hmax, abstol, reltol) printmessages && println("hinit = $h") nfcn += 1 end @@ -152,7 +152,7 @@ function dop853(F!::Function, p::ODEProblem, y0, tspan; F!(p, k1, x, y) end - y1, err = dopcore(n, F!, p, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, abstol, reltol) + y1, err = dopcore(n, p, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, abstol, reltol) xph = x+h nfcn += 11 @@ -246,7 +246,7 @@ function denseout(ind, t, told, h, coeff, dense) end end -function hinit(n::Int64, F!::Function, p::ODEProblem, x::Float64, y::Vector, xend::Float64, posneg::Float64, f0::AbstractVector, f1::AbstractVector, y0::AbstractVector, iord::Int64, hmax::Float64, abstol::Vector{Float64}, reltol::Vector{Float64}) +function hinit(n::Int64, p::ODEProblem, x::Float64, y::Vector, xend::Float64, posneg::Float64, f0::AbstractVector, f1::AbstractVector, y0::AbstractVector, iord::Int64, hmax::Float64, abstol::Vector{Float64}, reltol::Vector{Float64}) dnf = 0.0 dny = 0.0 for i = 1:n @@ -279,7 +279,7 @@ function hinit(n::Int64, F!::Function, p::ODEProblem, x::Float64, y::Vector, xen return h*posneg end -function dopcore(n::Int64, F!::Function, p::ODEProblem, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector, reltol::Vector) +function dopcore(n::Int64, p::ODEProblem, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector, reltol::Vector) a21 = 5.26001519587677318785587544488e-2 a31 = 1.97250569845378994544595329183e-2 a32 = 5.91751709536136983633785987549e-2 From b893824eea2973237f4c44b794a755c8140a5b69 Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Wed, 25 Jun 2014 17:39:18 +0200 Subject: [PATCH 11/20] Fix tests. --- src/dop853.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dop853.jl b/src/dop853.jl index 77d905707..485ba0c7a 100644 --- a/src/dop853.jl +++ b/src/dop853.jl @@ -9,7 +9,7 @@ F!(p::ODEProblemFunction, y, x, t) = copy!(y, [p.f(x, t)]) dop853(f::Function, y0::Vector, tspan; args...) = dop853(ODEProblemFunction(f), y0, tspan; args...) function dop853(f::Function, y0::Number, tspan; args...) - tout, yout = dop853(F!, ODEProblemFunction(f), [y0], tspan; args...) + tout, yout = dop853(ODEProblemFunction(f), [y0], tspan; args...) return tout, vcat(yout...) end From f752738a37e6167c066fbcf0e54fdc631a5798df Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Thu, 26 Jun 2014 14:11:41 +0200 Subject: [PATCH 12/20] Cleanup. --- src/dop853.jl | 70 +++++++++++++++++++++++++-------------------------- 1 file changed, 35 insertions(+), 35 deletions(-) diff --git a/src/dop853.jl b/src/dop853.jl index 485ba0c7a..39be8b9e5 100644 --- a/src/dop853.jl +++ b/src/dop853.jl @@ -73,6 +73,30 @@ function dop853(p::ODEProblem, y0, tspan; d714 = 0.96324553959188282948394950600e+02 d715 = -0.39177261675615439165231486172e+02 d716 = -0.14972683625798562581422125276e+03 + a141 = 5.61675022830479523392909219681e-2 + a147 = 2.53500210216624811088794765333e-1 + a148 = -2.46239037470802489917441475441e-1 + a149 = -1.24191423263816360469010140626e-1 + a1410 = 1.5329179827876569731206322685e-1 + a1411 = 8.20105229563468988491666602057e-3 + a1412 = 7.56789766054569976138603589584e-3 + a1413 = -8.298e-3 + a151 = 3.18346481635021405060768473261e-2 + a156 = 2.83009096723667755288322961402e-2 + a157 = 5.35419883074385676223797384372e-2 + a158 = -5.49237485713909884646569340306e-2 + a1511 = -1.08347328697249322858509316994e-4 + a1512 = 3.82571090835658412954920192323e-4 + a1513 = -3.40465008687404560802977114492e-4 + a1514 = 1.41312443674632500278074618366e-1 + a161 = -4.28896301583791923408573538692e-1 + a166 = -4.69762141536116384314449447206e0 + a167 = 7.68342119606259904184240953878e0 + a168 = 4.06898981839711007970213554331e0 + a169 = 3.56727187455281109270669543021e-1 + a1613 = -1.39902416515901462129418009734e-3 + a1614 = 2.9475147891527723389556272149e0 + a1615 = -9.15095847217987001081870187138e0 irtrn = 0 @@ -92,17 +116,17 @@ function dop853(p::ODEProblem, y0, tspan; nstep = 0 naccpt = 0 nrejct = 0 - k1 = 0. * y0 - k2 = 0. * y0 - k3 = 0. * y0 - k4 = 0. * y0 - k5 = 0. * y0 - k6 = 0. * y0 - k7 = 0. * y0 - k8 = 0. * y0 - k9 = 0. * y0 - k10 = 0. * y0 - y = 0. * y0 + k1 = 0.*y0 + k2 = 0.*y0 + k3 = 0.*y0 + k4 = 0.*y0 + k5 = 0.*y0 + k6 = 0.*y0 + k7 = 0.*y0 + k8 = 0.*y0 + k9 = 0.*y0 + k10 = 0.*y0 + y = 0.*y0 copy!(y, y0) facold = 1e-4 @@ -330,30 +354,6 @@ function dopcore(n::Int64, p::ODEProblem, x::Float64, y::Vector, h::Float64, k1: a129 = -8.87285693353062954433549289258e0 a1210 = 1.23605671757943030647266201528e1 a1211 = 6.43392746015763530355970484046e-1 - a141 = 5.61675022830479523392909219681e-2 - a147 = 2.53500210216624811088794765333e-1 - a148 = -2.46239037470802489917441475441e-1 - a149 = -1.24191423263816360469010140626e-1 - a1410 = 1.5329179827876569731206322685e-1 - a1411 = 8.20105229563468988491666602057e-3 - a1412 = 7.56789766054569976138603589584e-3 - a1413 = -8.298e-3 - a151 = 3.18346481635021405060768473261e-2 - a156 = 2.83009096723667755288322961402e-2 - a157 = 5.35419883074385676223797384372e-2 - a158 = -5.49237485713909884646569340306e-2 - a1511 = -1.08347328697249322858509316994e-4 - a1512 = 3.82571090835658412954920192323e-4 - a1513 = -3.40465008687404560802977114492e-4 - a1514 = 1.41312443674632500278074618366e-1 - a161 = -4.28896301583791923408573538692e-1 - a166 = -4.69762141536116384314449447206e0 - a167 = 7.68342119606259904184240953878e0 - a168 = 4.06898981839711007970213554331e0 - a169 = 3.56727187455281109270669543021e-1 - a1613 = -1.39902416515901462129418009734e-3 - a1614 = 2.9475147891527723389556272149e0 - a1615 = -9.15095847217987001081870187138e0 b1 = 5.42937341165687622380535766363e-2 b6 = 4.45031289275240888144113950566e0 b7 = 1.89151789931450038304281599044e0 From a39aa73e860a7ba2c2e1788e7f5cfdd49e9c8f50 Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Fri, 27 Jun 2014 11:18:11 +0200 Subject: [PATCH 13/20] Add coeff type in preparation for DOPRI5. --- src/ODE.jl | 2 +- src/dop853.jl | 457 --------------------------------- src/odedop.jl | 693 ++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 694 insertions(+), 458 deletions(-) delete mode 100644 src/dop853.jl create mode 100644 src/odedop.jl diff --git a/src/ODE.jl b/src/ODE.jl index ee4cf9276..e185ab5e0 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -13,7 +13,7 @@ export ode23, ode4, ode45, ode4s, ode4ms, ode78, dop853, ODEProblem # oderosenbrock, ode4s, ode4s_kr, ode4s_s, # ode4ms, ode_ms -include("dop853.jl") +include("odedop.jl") #ODE23 Solve non-stiff differential equations. # diff --git a/src/dop853.jl b/src/dop853.jl deleted file mode 100644 index 39be8b9e5..000000000 --- a/src/dop853.jl +++ /dev/null @@ -1,457 +0,0 @@ -abstract ODEProblem - -type ODEProblemFunction <: ODEProblem - f::Function -end - -F!(p::ODEProblemFunction, y, x, t) = copy!(y, [p.f(x, t)]) - -dop853(f::Function, y0::Vector, tspan; args...) = dop853(ODEProblemFunction(f), y0, tspan; args...) - -function dop853(f::Function, y0::Number, tspan; args...) - tout, yout = dop853(ODEProblemFunction(f), [y0], tspan; args...) - return tout, vcat(yout...) -end - -function dop853(p::ODEProblem, y0, tspan; - reltol=[1e-6], abstol=[sqrt(eps())], - safe=0.9, fac1=0.333, fac2=6.0, beta=0.0, - maxstep=tspan[end]-tspan[1], initstep=0.0, - maxsteps=100000, printmessages=false, nstiff=1000, - iout=0, solout=s(x...)=return, dense=[1:length(y0)], - points=:all) - - - c14 = 0.1e+00 - c15 = 0.2e+00 - c16 = 0.777777777777777777777777777778e+00 - d41 = -0.84289382761090128651353491142e+01 - d46 = 0.56671495351937776962531783590e+00 - d47 = -0.30689499459498916912797304727e+01 - d48 = 0.23846676565120698287728149680e+01 - d49 = 0.21170345824450282767155149946e+01 - d410 = -0.87139158377797299206789907490e+00 - d411 = 0.22404374302607882758541771650e+01 - d412 = 0.63157877876946881815570249290e+00 - d413 = -0.88990336451333310820698117400e-01 - d414 = 0.18148505520854727256656404962e+02 - d415 = -0.91946323924783554000451984436e+01 - d416 = -0.44360363875948939664310572000e+01 - d51 = 0.10427508642579134603413151009e+02 - d56 = 0.24228349177525818288430175319e+03 - d57 = 0.16520045171727028198505394887e+03 - d58 = -0.37454675472269020279518312152e+03 - d59 = -0.22113666853125306036270938578e+02 - d510 = 0.77334326684722638389603898808e+01 - d511 = -0.30674084731089398182061213626e+02 - d512 = -0.93321305264302278729567221706e+01 - d513 = 0.15697238121770843886131091075e+02 - d514 = -0.31139403219565177677282850411e+02 - d515 = -0.93529243588444783865713862664e+01 - d516 = 0.35816841486394083752465898540e+02 - d61 = 0.19985053242002433820987653617e+02 - d66 = -0.38703730874935176555105901742e+03 - d67 = -0.18917813819516756882830838328e+03 - d68 = 0.52780815920542364900561016686e+03 - d69 = -0.11573902539959630126141871134e+02 - d610 = 0.68812326946963000169666922661e+01 - d611 = -0.10006050966910838403183860980e+01 - d612 = 0.77771377980534432092869265740e+00 - d613 = -0.27782057523535084065932004339e+01 - d614 = -0.60196695231264120758267380846e+02 - d615 = 0.84320405506677161018159903784e+02 - d616 = 0.11992291136182789328035130030e+02 - d71 = -0.25693933462703749003312586129e+02 - d76 = -0.15418974869023643374053993627e+03 - d77 = -0.23152937917604549567536039109e+03 - d78 = 0.35763911791061412378285349910e+03 - d79 = 0.93405324183624310003907691704e+02 - d710 = -0.37458323136451633156875139351e+02 - d711 = 0.10409964950896230045147246184e+03 - d712 = 0.29840293426660503123344363579e+02 - d713 = -0.43533456590011143754432175058e+02 - d714 = 0.96324553959188282948394950600e+02 - d715 = -0.39177261675615439165231486172e+02 - d716 = -0.14972683625798562581422125276e+03 - a141 = 5.61675022830479523392909219681e-2 - a147 = 2.53500210216624811088794765333e-1 - a148 = -2.46239037470802489917441475441e-1 - a149 = -1.24191423263816360469010140626e-1 - a1410 = 1.5329179827876569731206322685e-1 - a1411 = 8.20105229563468988491666602057e-3 - a1412 = 7.56789766054569976138603589584e-3 - a1413 = -8.298e-3 - a151 = 3.18346481635021405060768473261e-2 - a156 = 2.83009096723667755288322961402e-2 - a157 = 5.35419883074385676223797384372e-2 - a158 = -5.49237485713909884646569340306e-2 - a1511 = -1.08347328697249322858509316994e-4 - a1512 = 3.82571090835658412954920192323e-4 - a1513 = -3.40465008687404560802977114492e-4 - a1514 = 1.41312443674632500278074618366e-1 - a161 = -4.28896301583791923408573538692e-1 - a166 = -4.69762141536116384314449447206e0 - a167 = 7.68342119606259904184240953878e0 - a168 = 4.06898981839711007970213554331e0 - a169 = 3.56727187455281109270669543021e-1 - a1613 = -1.39902416515901462129418009734e-3 - a1614 = 2.9475147891527723389556272149e0 - a1615 = -9.15095847217987001081870187138e0 - - irtrn = 0 - - x = tspan[1] - xend = tspan[end] - tout = Array(typeof(tspan[1]),0) - yout = Array(typeof(y0*one(x)),0) - - if points == :all - push!(tout, x) - push!(yout, y0) - end - - n = length(y0) - nrd = length(dense) - nfcn = 0 - nstep = 0 - naccpt = 0 - nrejct = 0 - k1 = 0.*y0 - k2 = 0.*y0 - k3 = 0.*y0 - k4 = 0.*y0 - k5 = 0.*y0 - k6 = 0.*y0 - k7 = 0.*y0 - k8 = 0.*y0 - k9 = 0.*y0 - k10 = 0.*y0 - y = 0.*y0 - - copy!(y, y0) - facold = 1e-4 - expo1 = 1.0/8.0 - beta*0.2 - facc1 = 1.0/fac1 - facc2 = 1.0/fac2 - posneg = sign(xend-x) - abstol = length(abstol) == 1 ? fill(abstol[1], n) : abstol - reltol = length(reltol) == 1 ? fill(reltol[1], n) : reltol - last = false - hlamb = 0.0 - iasti = 0 - F!(p, k1, x, y) - nfcn += 1 - hmax = abs(maxstep) - nmax = maxsteps - h = initstep - iord = 8 - if h == 0.0 - h = hinit(n, p, x, y, xend, posneg, k1, k2, k3, iord, hmax, abstol, reltol) - printmessages && println("hinit = $h") - nfcn += 1 - end - reject = false - xold = x - if iout != 0 - #hout = 1.0 - irtrn = solout(naccpt+1, xold, x, y, con, icomp) - if irtrn < 0 - return y - end - end - - while true - if nstep > nmax - error("Exit at x=$x. More than nmax=$nmax steps needed.") - end - if 0.1*abs(h) <= abs(x)*eps() - error("Exit at x=$x. Step size too small, h=$h.") - end - if (x+1.01*h-xend)*posneg > 0.0 - h = xend-x - last = true - end - nstep += 1 - if irtrn >= 2 - F!(p, k1, x, y) - end - - y1, err = dopcore(n, p, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, abstol, reltol) - xph = x+h - nfcn += 11 - - fac11 = err^expo1 - fac = fac11/facold^beta - fac = max(facc2, min(facc1, fac/safe)) - hnew = h/fac - if err <= 1.0 - facold = max(err, 1e-4) - naccpt += 1 - F!(p, k4, xph, k5) - nfcn += 1 - # Stiffness detection - if mod(naccpt, nstiff) == 0 || iasti > 0 - nonsti = 0 - stnum = 0.0 - stden = 0.0 - for i = 1:n - stnum += (k4[i] - k3[i])^2 - stden += (k5[i] - y1[i])^2 - end - if stden > 0.0 - hlamb = abs(h)*sqrt(stnum/stden) - end - if hlamb > 6.1 - nonsti = 0 - iasti += 1 - if iasti == 15 - if printmessages - println("The problem seems to become stiff at x=$x") - end - end - else - nonsti += 1 - if nonsti == 6 - iasti = 0 - end - end - end - # Final preparation for dense output - event = iout == 3 && xout <= xph - if iout == 2 || event - for (i,j) in enumerate(dense) - end - end - copy!(k1, k4) - copy!(y, k5) - xold = x - x = xph - if points == :all - push!(tout, x) - push!(yout, y1) - end - # if - # solout - # end - # Normal exit - if last - h = hnew - idid = 1 - if points == :last - return x, y1 - end - return tout, yout - end - if abs(hnew) > hmax - hnew = posneg*hmax - end - if reject - hnew = posneg*min(abs(hnew),abs(h)) - end - reject = false - else - hnew = h/min(facc1,fac11/safe) - reject = true - if naccpt >= 1 - nrejct += 1 - end - last = false - end - h = hnew - end -end - -function denseout(ind, t, told, h, coeff, dense) - if ~in(ind, dense) - error("No dense output available for component: $ind.") - else - s = (t - told)/h - s1 = 1.0 - s - end -end - -function hinit(n::Int64, p::ODEProblem, x::Float64, y::Vector, xend::Float64, posneg::Float64, f0::AbstractVector, f1::AbstractVector, y0::AbstractVector, iord::Int64, hmax::Float64, abstol::Vector{Float64}, reltol::Vector{Float64}) - dnf = 0.0 - dny = 0.0 - for i = 1:n - sk = abstol[i] + reltol[i]*abs(y[i]) - dnf += (f0[i]/sk)^2 - dny += (y[i]/sk)^2 - end - if maximum(abs(dnf)) <= 1e-10 || maximum(abs(dny)) <= 1e-10 - h = 1e-6 - else - h = sqrt(maximum(abs(dny/dnf)))*0.01 - end - h = min(h, hmax) - h = h*posneg - y1 = y + h*y0 - F!(p, f1, x+h, y1) - der2 = 0.0 - for i = 1:n - sk = abstol[i] + reltol[i]*abs(y[i]) - der2 += ((f1[i]-f0[i])/sk)^2 - end - der2 = sqrt(maximum(abs(der2)))/h - der12 = max(maximum(abs(der2)), sqrt(maximum(abs(dnf)))) - if der12 <= 1e-15 - h1 = max(1e-6, abs(h)*1e-3) - else - h1 = (0.01/der12)^(1.0/iord) - end - h = min(100*abs(h), h1, hmax) - return h*posneg -end - -function dopcore(n::Int64, p::ODEProblem, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector, reltol::Vector) - a21 = 5.26001519587677318785587544488e-2 - a31 = 1.97250569845378994544595329183e-2 - a32 = 5.91751709536136983633785987549e-2 - a41 = 2.95875854768068491816892993775e-2 - a43 = 8.87627564304205475450678981324e-2 - a51 = 2.41365134159266685502369798665e-1 - a53 = -8.84549479328286085344864962717e-1 - a54 = 9.24834003261792003115737966543e-1 - a61 = 3.7037037037037037037037037037e-2 - a64 = 1.70828608729473871279604482173e-1 - a65 = 1.25467687566822425016691814123e-1 - a71 = 3.7109375e-2 - a74 = 1.70252211019544039314978060272e-1 - a75 = 6.02165389804559606850219397283e-2 - a76 = -1.7578125e-2 - a81 = 3.70920001185047927108779319836e-2 - a84 = 1.70383925712239993810214054705e-1 - a85 = 1.07262030446373284651809199168e-1 - a86 = -1.53194377486244017527936158236e-2 - a87 = 8.27378916381402288758473766002e-3 - a91 = 6.24110958716075717114429577812e-1 - a94 = -3.36089262944694129406857109825e0 - a95 = -8.68219346841726006818189891453e-1 - a96 = 2.75920996994467083049415600797e1 - a97 = 2.01540675504778934086186788979e1 - a98 = -4.34898841810699588477366255144e1 - a101 = 4.77662536438264365890433908527e-1 - a104 = -2.48811461997166764192642586468e0 - a105 = -5.90290826836842996371446475743e-1 - a106 = 2.12300514481811942347288949897e1 - a107 = 1.52792336328824235832596922938e1 - a108 = -3.32882109689848629194453265587e1 - a109 = -2.03312017085086261358222928593e-2 - a111 = -9.3714243008598732571704021658e-1 - a114 = 5.18637242884406370830023853209e0 - a115 = 1.09143734899672957818500254654e0 - a116 = -8.14978701074692612513997267357e0 - a117 = -1.85200656599969598641566180701e1 - a118 = 2.27394870993505042818970056734e1 - a119 = 2.49360555267965238987089396762e0 - a1110 = -3.0467644718982195003823669022e0 - a121 = 2.27331014751653820792359768449e0 - a124 = -1.05344954667372501984066689879e1 - a125 = -2.00087205822486249909675718444e0 - a126 = -1.79589318631187989172765950534e1 - a127 = 2.79488845294199600508499808837e1 - a128 = -2.85899827713502369474065508674e0 - a129 = -8.87285693353062954433549289258e0 - a1210 = 1.23605671757943030647266201528e1 - a1211 = 6.43392746015763530355970484046e-1 - b1 = 5.42937341165687622380535766363e-2 - b6 = 4.45031289275240888144113950566e0 - b7 = 1.89151789931450038304281599044e0 - b8 = -5.8012039600105847814672114227e0 - b9 = 3.1116436695781989440891606237e-1 - b10 = -1.52160949662516078556178806805e-1 - b11 = 2.01365400804030348374776537501e-1 - b12 = 4.47106157277725905176885569043e-2 - c2 = 0.526001519587677318785587544488e-01 - c3 = 0.789002279381515978178381316732e-01 - c4 = 0.118350341907227396726757197510e+00 - c5 = 0.281649658092772603273242802490e+00 - c6 = 0.333333333333333333333333333333e+00 - c7 = 0.25e+00 - c8 = 0.307692307692307692307692307692e+00 - c9 = 0.651282051282051282051282051282e+00 - c10 = 0.6e+00 - c11 = 0.857142857142857142857142857142e+00 - bhh1 = 0.244094488188976377952755905512e+00 - bhh2 = 0.733846688281611857341361741547e+00 - bhh3 = 0.220588235294117647058823529412e-01 - er1 = 0.1312004499419488073250102996e-01 - er6 = -0.1225156446376204440720569753e+01 - er7 = -0.4957589496572501915214079952e+00 - er8 = 0.1664377182454986536961530415e+01 - er9 = -0.3503288487499736816886487290e+00 - er10 = 0.3341791187130174790297318841e+00 - er11 = 0.8192320648511571246570742613e-01 - er12 = -0.2235530786388629525884427845e-01 - - y1 = 0.0 * y - - for i = 1:n - y1[i] = y[i] + h*a21*k1[i] - end - F!(p, k2, x+c2*h, y1) - for i = 1:n - y1[i] = y[i] + h*(a31*k1[i] + a32*k2[i]) - end - F!(p, k3, x+c3*h, y1) - for i = 1:n - y1[i] = y[i]+h*(a41*k1[i]+a43*k3[i]) - end - F!(p, k4, x+c4*h, y1) - for i = 1:n - y1[i] = y[i]+h*(a51*k1[i]+a53*k3[i]+a54*k4[i]) - end - F!(p, k5, x+c5*h, y1) - for i = 1:n - y1[i] = y[i]+h*(a61*k1[i]+a64*k4[i]+a65*k5[i]) - end - F!(p, k6, x+c6*h, y1) - for i = 1:n - y1[i] = y[i]+h*(a71*k1[i]+a74*k4[i]+a75*k5[i]+a76*k6[i]) - end - F!(p, k7, x+c7*h, y1) - for i = 1:n - y1[i] = y[i]+h*(a81*k1[i]+a84*k4[i]+a85*k5[i]+a86*k6[i]+a87*k7[i]) - end - F!(p, k8, x+c8*h, y1) - for i = 1:n - y1[i] = y[i]+h*(a91*k1[i]+a94*k4[i]+a95*k5[i]+a96*k6[i]+a97*k7[i]+a98*k8[i]) - end - F!(p, k9, x+c9*h, y1) - for i = 1:n - y1[i] = y[i]+h*(a101*k1[i]+a104*k4[i]+a105*k5[i]+a106*k6[i] - +a107*k7[i]+a108*k8[i]+a109*k9[i]) - end - F!(p, k10, x+c10*h, y1) - for i = 1:n - y1[i] = y[i]+h*(a111*k1[i]+a114*k4[i]+a115*k5[i]+a116*k6[i] - +a117*k7[i]+a118*k8[i]+a119*k9[i]+a1110*k10[i]) - end - F!(p, k2, x+c11*h, y1) - for i = 1:n - y1[i] = y[i]+h*(a121*k1[i]+a124*k4[i]+a125*k5[i]+a126*k6[i] - +a127*k7[i]+a128*k8[i]+a129*k9[i]+a1210*k10[i]+a1211*k2[i]) - end - F!(p, k3, x+h, y1) - for i = 1:n - k4[i] = b1*k1[i]+b6*k6[i]+b7*k7[i]+b8*k8[i]+b9*k9[i]+b10*k10[i]+b11*k2[i]+b12*k3[i] - k5[i] = y[i]+h*k4[i] - end - - # Error estimation - err = 0.0 - err2 = 0.0 - for i = 1:n - sk = abstol[i] + reltol[i]*max(maximum(abs(y[i])), maximum(abs(k5[i]))) - erri = k4[i] - bhh1*k1[i] - bhh2*k9[i] - bhh3*k3[i] - err2 += (erri/sk)*(erri/sk) - erri = er1*k1[i] + er6*k6[i] + er7*k7[i] + er8*k8[i] + er9*k9[i] + er10*k10[i] + er11*k2[i] + er12*k3[i] - err += (erri/sk)*(erri/sk) - end - deno = maximum(abs(err)) + 0.01*maximum(abs(err2)) - if deno <= 0.0 - deno = 1.0 - end - err = abs(h)*err*sqrt(1.0/(n*deno)) - return y1, maximum(abs(err)) -end diff --git a/src/odedop.jl b/src/odedop.jl new file mode 100644 index 000000000..3be1d1c84 --- /dev/null +++ b/src/odedop.jl @@ -0,0 +1,693 @@ +abstract ODEProblem + +type ODEProblemFunction + f::Function +end + +type DOP853 + a21::Float64 + a31::Float64 + a32::Float64 + a41::Float64 + a43::Float64 + a51::Float64 + a53::Float64 + a54::Float64 + a61::Float64 + a64::Float64 + a65::Float64 + a71::Float64 + a74::Float64 + a75::Float64 + a76::Float64 + a81::Float64 + a84::Float64 + a85::Float64 + a86::Float64 + a87::Float64 + a91::Float64 + a94::Float64 + a95::Float64 + a96::Float64 + a97::Float64 + a98::Float64 + a101::Float64 + a104::Float64 + a105::Float64 + a106::Float64 + a107::Float64 + a108::Float64 + a109::Float64 + a111::Float64 + a114::Float64 + a115::Float64 + a116::Float64 + a117::Float64 + a118::Float64 + a119::Float64 + a1110::Float64 + a121::Float64 + a124::Float64 + a125::Float64 + a126::Float64 + a127::Float64 + a128::Float64 + a129::Float64 + a1210::Float64 + a1211::Float64 + b1::Float64 + b6::Float64 + b7::Float64 + b8::Float64 + b9::Float64 + b10::Float64 + b11::Float64 + b12::Float64 + c2::Float64 + c3::Float64 + c4::Float64 + c5::Float64 + c6::Float64 + c7::Float64 + c8::Float64 + c9::Float64 + c10::Float64 + c11::Float64 + bhh1::Float64 + bhh2::Float64 + bhh3::Float64 + er1::Float64 + er6::Float64 + er7::Float64 + er8::Float64 + er9::Float64 + er10::Float64 + er11::Float64 + er12::Float64 + c14::Float64 + c15::Float64 + c16::Float64 + d41::Float64 + d46::Float64 + d47::Float64 + d48::Float64 + d49::Float64 + d410::Float64 + d411::Float64 + d412::Float64 + d413::Float64 + d414::Float64 + d415::Float64 + d416::Float64 + d51::Float64 + d56::Float64 + d57::Float64 + d58::Float64 + d59::Float64 + d510::Float64 + d511::Float64 + d512::Float64 + d513::Float64 + d514::Float64 + d515::Float64 + d516::Float64 + d61::Float64 + d66::Float64 + d67::Float64 + d68::Float64 + d69::Float64 + d610::Float64 + d611::Float64 + d612::Float64 + d613::Float64 + d614::Float64 + d615::Float64 + d616::Float64 + d71::Float64 + d76::Float64 + d77::Float64 + d78::Float64 + d79::Float64 + d710::Float64 + d711::Float64 + d712::Float64 + d713::Float64 + d714::Float64 + d715::Float64 + d716::Float64 + a141::Float64 + a147::Float64 + a148::Float64 + a149::Float64 + a1410::Float64 + a1411::Float64 + a1412::Float64 + a1413::Float64 + a151::Float64 + a156::Float64 + a157::Float64 + a158::Float64 + a1511::Float64 + a1512::Float64 + a1513::Float64 + a1514::Float64 + a161::Float64 + a166::Float64 + a167::Float64 + a168::Float64 + a169::Float64 + a1613::Float64 + a1614::Float64 + a1615::Float64 +end + +type DOPRI5 + c2::Float64 + c3::Float64 + c4::Float64 + c5::Float64 + a21::Float64 + a31::Float64 + a32::Float64 + a41::Float64 + a42::Float64 + a43::Float64 + a51::Float64 + a52::Float64 + a53::Float64 + a54::Float64 + a61::Float64 + a62::Float64 + a63::Float64 + a64::Float64 + a65::Float64 + a71::Float64 + a73::Float64 + a74::Float64 + a75::Float64 + a76::Float64 + e1::Float64 + e3::Float64 + e4::Float64 + e5::Float64 + e6::Float64 + e7::Float64 + d1::Float64 + d3::Float64 + d4::Float64 + d5::Float64 + d6::Float64 + d7::Float64 +end + +const dopri5coeff = DOPRI5( + 0.2, + 0.3, + 0.8, + 8.0/9.0, + 0.2, + 3.0/40.0, + 9.0/40.0, + 44.0/45.0, + -56.0/15.0, + 32.0/9.0, + 19372.0/6561.0, + -25360.0/2187.0, + 64448.0/6561.0, + -212.0/729.0, + 9017.0/3168.0, + -355.0/33.0, + 46732.0/5247.0, + 49.0/176.0, + -5103.0/18656.0, + 35.0/384.0, + 500.0/1113.0, + 125.0/192.0, + -2187.0/6784.0, + 11.0/84.0, + 71.0/57600.0, + -71.0/16695.0, + 71.0/1920.0, + -17253.0/339200.0, + 22.0/525.0, + -1.0/40.0 , + -12715105075.0/11282082432.0, + 87487479700.0/32700410799.0, + -10690763975.0/1880347072.0, + 701980252875.0/199316789632.0, + -1453857185.0/822651844.0, + 69997945.0/29380423.0) + +const dop853coeff = DOP853(5.26001519587677318785587544488e-2, + 1.97250569845378994544595329183e-2, + 5.91751709536136983633785987549e-2, + 2.95875854768068491816892993775e-2, + 8.87627564304205475450678981324e-2, + 2.41365134159266685502369798665e-1, + -8.84549479328286085344864962717e-1, + 9.24834003261792003115737966543e-1, + 3.7037037037037037037037037037e-2, + 1.70828608729473871279604482173e-1, + 1.25467687566822425016691814123e-1, + 3.7109375e-2, + 1.70252211019544039314978060272e-1, + 6.02165389804559606850219397283e-2, + -1.7578125e-2, + 3.70920001185047927108779319836e-2, + 1.70383925712239993810214054705e-1, + 1.07262030446373284651809199168e-1, + -1.53194377486244017527936158236e-2, + 8.27378916381402288758473766002e-3, + 6.24110958716075717114429577812e-1, + -3.36089262944694129406857109825e0, + -8.68219346841726006818189891453e-1, + 2.75920996994467083049415600797e1, + 2.01540675504778934086186788979e1, + -4.34898841810699588477366255144e1, + 4.77662536438264365890433908527e-1, + -2.48811461997166764192642586468e0, + -5.90290826836842996371446475743e-1, + 2.12300514481811942347288949897e1, + 1.52792336328824235832596922938e1, + -3.32882109689848629194453265587e1, + -2.03312017085086261358222928593e-2, + -9.3714243008598732571704021658e-1, + 5.18637242884406370830023853209e0, + 1.09143734899672957818500254654e0, + -8.14978701074692612513997267357e0, + -1.85200656599969598641566180701e1, + 2.27394870993505042818970056734e1, + 2.49360555267965238987089396762e0, + -3.0467644718982195003823669022e0, + 2.27331014751653820792359768449e0, + -1.05344954667372501984066689879e1, + -2.00087205822486249909675718444e0, + -1.79589318631187989172765950534e1, + 2.79488845294199600508499808837e1, + -2.85899827713502369474065508674e0, + -8.87285693353062954433549289258e0, + 1.23605671757943030647266201528e1, + 6.43392746015763530355970484046e-1, + 5.42937341165687622380535766363e-2, + 4.45031289275240888144113950566e0, + 1.89151789931450038304281599044e0, + -5.8012039600105847814672114227e0, + 3.1116436695781989440891606237e-1, + -1.52160949662516078556178806805e-1, + 2.01365400804030348374776537501e-1, + 4.47106157277725905176885569043e-2, + 0.526001519587677318785587544488e-01, + 0.789002279381515978178381316732e-01, + 0.118350341907227396726757197510e+00, + 0.281649658092772603273242802490e+00, + 0.333333333333333333333333333333e+00, + 0.25e+00, + 0.307692307692307692307692307692e+00, + 0.651282051282051282051282051282e+00, + 0.6e+00, + 0.857142857142857142857142857142e+00, + 0.244094488188976377952755905512e+00, + 0.733846688281611857341361741547e+00, + 0.220588235294117647058823529412e-01, + 0.1312004499419488073250102996e-01, + -0.1225156446376204440720569753e+01, + -0.4957589496572501915214079952e+00, + 0.1664377182454986536961530415e+01, + -0.3503288487499736816886487290e+00, + 0.3341791187130174790297318841e+00, + 0.8192320648511571246570742613e-01, + -0.2235530786388629525884427845e-01, + 0.1e+00, + 0.2e+00, + 0.777777777777777777777777777778e+00, + -0.84289382761090128651353491142e+01, + 0.56671495351937776962531783590e+00, + -0.30689499459498916912797304727e+01, + 0.23846676565120698287728149680e+01, + 0.21170345824450282767155149946e+01, + -0.87139158377797299206789907490e+00, + 0.22404374302607882758541771650e+01, + 0.63157877876946881815570249290e+00, + -0.88990336451333310820698117400e-01, + 0.18148505520854727256656404962e+02, + -0.91946323924783554000451984436e+01, + -0.44360363875948939664310572000e+01, + 0.10427508642579134603413151009e+02, + 0.24228349177525818288430175319e+03, + 0.16520045171727028198505394887e+03, + -0.37454675472269020279518312152e+03, + -0.22113666853125306036270938578e+02, + 0.77334326684722638389603898808e+01, + -0.30674084731089398182061213626e+02, + -0.93321305264302278729567221706e+01, + 0.15697238121770843886131091075e+02, + -0.31139403219565177677282850411e+02, + -0.93529243588444783865713862664e+01, + 0.35816841486394083752465898540e+02, + 0.19985053242002433820987653617e+02, + -0.38703730874935176555105901742e+03, + -0.18917813819516756882830838328e+03, + 0.52780815920542364900561016686e+03, + -0.11573902539959630126141871134e+02, + 0.68812326946963000169666922661e+01, + -0.10006050966910838403183860980e+01, + 0.77771377980534432092869265740e+00, + -0.27782057523535084065932004339e+01, + -0.60196695231264120758267380846e+02, + 0.84320405506677161018159903784e+02, + 0.11992291136182789328035130030e+02, + -0.25693933462703749003312586129e+02, + -0.15418974869023643374053993627e+03, + -0.23152937917604549567536039109e+03, + 0.35763911791061412378285349910e+03, + 0.93405324183624310003907691704e+02, + -0.37458323136451633156875139351e+02, + 0.10409964950896230045147246184e+03, + 0.29840293426660503123344363579e+02, + -0.43533456590011143754432175058e+02, + 0.96324553959188282948394950600e+02, + -0.39177261675615439165231486172e+02, + -0.14972683625798562581422125276e+03, + 5.61675022830479523392909219681e-2, + 2.53500210216624811088794765333e-1, + -2.46239037470802489917441475441e-1, + -1.24191423263816360469010140626e-1, + 1.5329179827876569731206322685e-1, + 8.20105229563468988491666602057e-3, + 7.56789766054569976138603589584e-3, + -8.298e-3, + 3.18346481635021405060768473261e-2, + 2.83009096723667755288322961402e-2, + 5.35419883074385676223797384372e-2, + -5.49237485713909884646569340306e-2, + -1.08347328697249322858509316994e-4, + 3.82571090835658412954920192323e-4, + -3.40465008687404560802977114492e-4, + 1.41312443674632500278074618366e-1, + -4.28896301583791923408573538692e-1, + -4.69762141536116384314449447206e0, + 7.68342119606259904184240953878e0, + 4.06898981839711007970213554331e0, + 3.56727187455281109270669543021e-1, + -1.39902416515901462129418009734e-3, + 2.9475147891527723389556272149e0, + -9.15095847217987001081870187138e0) + +F!(p::ODEProblemFunction, y, x, t) = copy!(y, [p.f(x, t)]) + +dop853(f::Function, y0::Vector, tspan; args...) = dop853(ODEProblemFunction(f), y0, tspan; args...) + +function dop853(f::Function, y0::Number, tspan; args...) + tout, yout = dop853(ODEProblemFunction(f), [y0], tspan; args...) + return tout, vcat(yout...) +end + + +dop853(p, y0, tspan; args...) = dop853(dop853coeff, p, y0, tspan; args...) + +function dop853(coeff, p, y0, tspan; + reltol=[1e-6], abstol=[sqrt(eps())], + safe=0.9, fac1=0.333, fac2=6.0, beta=0.0, + maxstep=tspan[end]-tspan[1], initstep=0.0, + maxsteps=100000, printmessages=false, nstiff=1000, + iout=0, solout=s(x...)=return, dense=[1:length(y0)], + points=:all) + + irtrn = 0 + + x = tspan[1] + xend = tspan[end] + tout = Array(typeof(tspan[1]),0) + yout = Array(typeof(y0*one(x)),0) + + if points == :all + push!(tout, x) + push!(yout, y0) + end + + n = length(y0) + nrd = length(dense) + nfcn = 0 + nstep = 0 + naccpt = 0 + nrejct = 0 + k1 = 0.*y0 + k2 = 0.*y0 + k3 = 0.*y0 + k4 = 0.*y0 + k5 = 0.*y0 + k6 = 0.*y0 + k7 = 0.*y0 + k8 = 0.*y0 + k9 = 0.*y0 + k10 = 0.*y0 + y = 0.*y0 + + copy!(y, y0) + facold = 1e-4 + expo1 = 1.0/8.0 - beta*0.2 + facc1 = 1.0/fac1 + facc2 = 1.0/fac2 + posneg = sign(xend-x) + abstol = length(abstol) == 1 ? fill(abstol[1], n) : abstol + reltol = length(reltol) == 1 ? fill(reltol[1], n) : reltol + last = false + hlamb = 0.0 + iasti = 0 + F!(p, k1, x, y) + nfcn += 1 + hmax = abs(maxstep) + nmax = maxsteps + h = initstep + iord = 8 + if h == 0.0 + h = hinit(n, p, x, y, xend, posneg, k1, k2, k3, iord, hmax, abstol, reltol) + printmessages && println("hinit = $h") + nfcn += 1 + end + reject = false + xold = x + if iout != 0 + #hout = 1.0 + irtrn = solout(naccpt+1, xold, x, y, con, icomp) + if irtrn < 0 + return y + end + end + + while true + if nstep > nmax + error("Exit at x=$x. More than nmax=$nmax steps needed.") + end + if 0.1*abs(h) <= abs(x)*eps() + error("Exit at x=$x. Step size too small, h=$h.") + end + if (x+1.01*h-xend)*posneg > 0.0 + h = xend-x + last = true + end + nstep += 1 + if irtrn >= 2 + F!(p, k1, x, y) + end + + y1, err = dopcore(coeff, n, p, x, y, h, k1, k2, k3, k4, k5, k6, k7, + k8, k9, k10, abstol, reltol) + xph = x+h + nfcn += 11 + + fac11 = err^expo1 + fac = fac11/facold^beta + fac = max(facc2, min(facc1, fac/safe)) + hnew = h/fac + if err <= 1.0 + facold = max(err, 1e-4) + naccpt += 1 + F!(p, k4, xph, k5) + nfcn += 1 + # Stiffness detection + if mod(naccpt, nstiff) == 0 || iasti > 0 + nonsti = 0 + stnum = 0.0 + stden = 0.0 + for i = 1:n + stnum += (k4[i] - k3[i])^2 + stden += (k5[i] - y1[i])^2 + end + if stden > 0.0 + hlamb = abs(h)*sqrt(stnum/stden) + end + if hlamb > 6.1 + nonsti = 0 + iasti += 1 + if iasti == 15 + if printmessages + println("The problem seems to become stiff at x=$x") + end + end + else + nonsti += 1 + if nonsti == 6 + iasti = 0 + end + end + end + # Final preparation for dense output + event = iout == 3 && xout <= xph + if iout == 2 || event + for (i,j) in enumerate(dense) + end + end + copy!(k1, k4) + copy!(y, k5) + xold = x + x = xph + if points == :all + push!(tout, x) + push!(yout, y1) + end + # if + # solout + # end + # Normal exit + if last + h = hnew + idid = 1 + if points == :last + return x, y1 + end + return tout, yout + end + if abs(hnew) > hmax + hnew = posneg*hmax + end + if reject + hnew = posneg*min(abs(hnew),abs(h)) + end + reject = false + else + hnew = h/min(facc1,fac11/safe) + reject = true + if naccpt >= 1 + nrejct += 1 + end + last = false + end + h = hnew + end +end + +function denseout(ind, t, told, h, coeff, dense) + if ~in(ind, dense) + error("No dense output available for component: $ind.") + else + s = (t - told)/h + s1 = 1.0 - s + end +end + +function hinit(n::Int64, p, x::Float64, y::Vector, xend::Float64, posneg::Float64, f0::AbstractVector, f1::AbstractVector, y0::AbstractVector, iord::Int64, hmax::Float64, abstol::Vector{Float64}, reltol::Vector{Float64}) + dnf = 0.0 + dny = 0.0 + for i = 1:n + sk = abstol[i] + reltol[i]*abs(y[i]) + dnf += (f0[i]/sk)^2 + dny += (y[i]/sk)^2 + end + if maximum(abs(dnf)) <= 1e-10 || maximum(abs(dny)) <= 1e-10 + h = 1e-6 + else + h = sqrt(maximum(abs(dny/dnf)))*0.01 + end + h = min(h, hmax) + h = h*posneg + y1 = y + h*y0 + F!(p, f1, x+h, y1) + der2 = 0.0 + for i = 1:n + sk = abstol[i] + reltol[i]*abs(y[i]) + der2 += ((f1[i]-f0[i])/sk)^2 + end + der2 = sqrt(maximum(abs(der2)))/h + der12 = max(maximum(abs(der2)), sqrt(maximum(abs(dnf)))) + if der12 <= 1e-15 + h1 = max(1e-6, abs(h)*1e-3) + else + h1 = (0.01/der12)^(1.0/iord) + end + h = min(100*abs(h), h1, hmax) + return h*posneg +end + +function dopcore(c::DOP853, n::Int64, p, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector, reltol::Vector) + y1 = 0.0 * y + + for i = 1:n + y1[i] = y[i] + h*c.a21*k1[i] + end + F!(p, k2, x+c.c2*h, y1) + for i = 1:n + y1[i] = y[i] + h*(c.a31*k1[i] + c.a32*k2[i]) + end + F!(p, k3, x+c.c3*h, y1) + for i = 1:n + y1[i] = y[i]+h*(c.a41*k1[i]+c.a43*k3[i]) + end + F!(p, k4, x+c.c4*h, y1) + for i = 1:n + y1[i] = y[i]+h*(c.a51*k1[i]+c.a53*k3[i]+c.a54*k4[i]) + end + F!(p, k5, x+c.c5*h, y1) + for i = 1:n + y1[i] = y[i]+h*(c.a61*k1[i]+c.a64*k4[i]+c.a65*k5[i]) + end + F!(p, k6, x+c.c6*h, y1) + for i = 1:n + y1[i] = y[i]+h*(c.a71*k1[i]+c.a74*k4[i]+c.a75*k5[i]+c.a76*k6[i]) + end + F!(p, k7, x+c.c7*h, y1) + for i = 1:n + y1[i] = y[i]+h*(c.a81*k1[i]+c.a84*k4[i]+c.a85*k5[i]+c.a86*k6[i]+c.a87*k7[i]) + end + F!(p, k8, x+c.c8*h, y1) + for i = 1:n + y1[i] = y[i]+h*(c.a91*k1[i]+c.a94*k4[i]+c.a95*k5[i]+c.a96*k6[i]+c.a97*k7[i]+c.a98*k8[i]) + end + F!(p, k9, x+c.c9*h, y1) + for i = 1:n + y1[i] = y[i]+h*(c.a101*k1[i]+c.a104*k4[i]+c.a105*k5[i]+c.a106*k6[i] + +c.a107*k7[i]+c.a108*k8[i]+c.a109*k9[i]) + end + F!(p, k10, x+c.c10*h, y1) + for i = 1:n + y1[i] = y[i]+h*(c.a111*k1[i]+c.a114*k4[i]+c.a115*k5[i]+c.a116*k6[i] + +c.a117*k7[i]+c.a118*k8[i]+c.a119*k9[i]+c.a1110*k10[i]) + end + F!(p, k2, x+c.c11*h, y1) + for i = 1:n + y1[i] = y[i]+h*(c.a121*k1[i]+c.a124*k4[i]+c.a125*k5[i]+c.a126*k6[i] + +c.a127*k7[i]+c.a128*k8[i]+c.a129*k9[i]+c.a1210*k10[i]+c.a1211*k2[i]) + end + F!(p, k3, x+h, y1) + for i = 1:n + k4[i] = c.b1*k1[i]+c.b6*k6[i]+c.b7*k7[i]+c.b8*k8[i]+c.b9*k9[i]+c.b10*k10[i]+c.b11*k2[i]+c.b12*k3[i] + k5[i] = y[i]+h*k4[i] + end + + # Error estimation + err = 0.0 + err2 = 0.0 + for i = 1:n + sk = abstol[i] + reltol[i]*max(maximum(abs(y[i])), maximum(abs(k5[i]))) + erri = k4[i] - c.bhh1*k1[i] - c.bhh2*k9[i] - c.bhh3*k3[i] + err2 += (erri/sk)*(erri/sk) + erri = c.er1*k1[i] + c.er6*k6[i] + c.er7*k7[i] + c.er8*k8[i] + c.er9*k9[i] + c.er10*k10[i] + c.er11*k2[i] + c.er12*k3[i] + err += (erri/sk)*(erri/sk) + end + deno = maximum(abs(err)) + 0.01*maximum(abs(err2)) + if deno <= 0.0 + deno = 1.0 + end + err = abs(h)*err*sqrt(1.0/(n*deno)) + return y1, maximum(abs(err)) +end From c4e703fba8c107d1050a785eb57bf94406ff029f Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Fri, 27 Jun 2014 13:21:37 +0200 Subject: [PATCH 14/20] Implement DOPRI5. --- src/ODE.jl | 2 +- src/odedop.jl | 83 ++++++++++++++++++++++++++++++++++++++++-------- test/runtests.jl | 3 +- 3 files changed, 73 insertions(+), 15 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index e185ab5e0..f5f1dcf00 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -5,7 +5,7 @@ module ODE using Polynomial ## minimal function export list -export ode23, ode4, ode45, ode4s, ode4ms, ode78, dop853, ODEProblem +export ode23, ode4, ode45, ode4s, ode4ms, ode78, dop853, dopri5, ODEProblem ## complete function export list #export ode23, ode4, diff --git a/src/odedop.jl b/src/odedop.jl index 3be1d1c84..ac1aa8286 100644 --- a/src/odedop.jl +++ b/src/odedop.jl @@ -1,9 +1,3 @@ -abstract ODEProblem - -type ODEProblemFunction - f::Function -end - type DOP853 a21::Float64 a31::Float64 @@ -393,19 +387,31 @@ const dop853coeff = DOP853(5.26001519587677318785587544488e-2, 2.9475147891527723389556272149e0, -9.15095847217987001081870187138e0) +abstract ODEProblem + +type ODEProblemFunction + f::Function +end + F!(p::ODEProblemFunction, y, x, t) = copy!(y, [p.f(x, t)]) +dop853(p, y0, tspan; args...) = odedop(dop853coeff, p, y0, tspan; args...) +dopri5(p, y0, tspan; args...) = odedop(dopri5coeff, p, y0, tspan; args...) + dop853(f::Function, y0::Vector, tspan; args...) = dop853(ODEProblemFunction(f), y0, tspan; args...) +dopri5(f::Function, y0::Vector, tspan; args...) = dopri5(ODEProblemFunction(f), y0, tspan; args...) function dop853(f::Function, y0::Number, tspan; args...) tout, yout = dop853(ODEProblemFunction(f), [y0], tspan; args...) return tout, vcat(yout...) end +function dopri5(f::Function, y0::Number, tspan; args...) + tout, yout = dopri5(ODEProblemFunction(f), [y0], tspan; args...) + return tout, vcat(yout...) +end -dop853(p, y0, tspan; args...) = dop853(dop853coeff, p, y0, tspan; args...) - -function dop853(coeff, p, y0, tspan; +function odedop(coeff, p, y0, tspan; reltol=[1e-6], abstol=[sqrt(eps())], safe=0.9, fac1=0.333, fac2=6.0, beta=0.0, maxstep=tspan[end]-tspan[1], initstep=0.0, @@ -511,8 +517,13 @@ function dop853(coeff, p, y0, tspan; stnum = 0.0 stden = 0.0 for i = 1:n - stnum += (k4[i] - k3[i])^2 - stden += (k5[i] - y1[i])^2 + if typeof(coeff) == DOP853 + stnum += (k4[i] - k3[i])^2 + stden += (k5[i] - y1[i])^2 + elseif typeof(coeff) == DOPRI5 + stnum += (k2[i] - k6[i])^2 + stden += (y1[i] - k7[i])^2 + end end if stden > 0.0 hlamb = abs(h)*sqrt(stnum/stden) @@ -538,8 +549,13 @@ function dop853(coeff, p, y0, tspan; for (i,j) in enumerate(dense) end end - copy!(k1, k4) - copy!(y, k5) + if typeof(coeff) == DOP853 + copy!(k1, k4) + copy!(y, k5) + elseif typeof(coeff) == DOPRI5 + copy!(k1, k2) + copy!(y, y1) + end xold = x x = xph if points == :all @@ -691,3 +707,44 @@ function dopcore(c::DOP853, n::Int64, p, x::Float64, y::Vector, h::Float64, k1:: err = abs(h)*err*sqrt(1.0/(n*deno)) return y1, maximum(abs(err)) end + +function dopcore(c::DOPRI5, n::Int64, p, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector, reltol::Vector) + y1 = 0.0 * y + + for i = 1:n + y1[i] = y[i] + h*c.a21*k1[i] + end + F!(p, k2, x+c.c2*h, y1) + for i = 1:n + y1[i] = y[i] + h*(c.a31*k1[i] + c.a32*k2[i]) + end + F!(p, k3, x+c.c3*h, y1) + for i = 1:n + y1[i] = y[i]+h*(c.a41*k1[i]+c.a42*k2[i]+c.a43*k3[i]) + end + F!(p, k4, x+c.c4*h, y1) + for i = 1:n + y1[i] = y[i]+h*(c.a51*k1[i]+c.a52*k2[i]+c.a53*k3[i]+c.a54*k4[i]) + end + F!(p, k5, x+c.c5*h, y1) + for i = 1:n + k7[i] = y[i]+h*(c.a61*k1[i]+c.a62*k2[i]+c.a63*k3[i]+c.a64*k4[i]+c.a65*k5[i]) + end + F!(p, k6, x+h, k7) + for i = 1:n + y1[i] = y[i]+h*(c.a71*k1[i]+c.a73*k3[i]+c.a74*k4[i]+c.a75*k5[i]+c.a76*k6[i]) + end + F!(p, k2, x+h, y1) + for i = 1:n + k4[i] = (c.e1*k1[i]+c.e3*k3[i]+c.e4*k4[i]+c.e5*k5[i]+c.e6*k6[i]+c.e7*k2[i])*h + end + + # Error estimation + err = 0.0 + for i = 1:n + sk = abstol[i] + reltol[i]*max(maximum(abs(y[i])), maximum(abs(y1[i]))) + err += (k4[i]/sk)*(k4[i]/sk) + end + err = sqrt(err/n) + return y1, maximum(abs(err)) +end diff --git a/test/runtests.jl b/test/runtests.jl index 2d55a401b..3ad28ae91 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,7 +17,8 @@ solvers = [ ODE.ode4s_kr, ODE.ode78_fb, - ODE.dop853] + ODE.dop853, + ODE.dopri5] for solver in solvers println("using $solver") From d8374ccfab83f953633c5735956a660aa5a0387b Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Fri, 27 Jun 2014 13:34:18 +0200 Subject: [PATCH 15/20] Minor adjustments. --- src/odedop.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/odedop.jl b/src/odedop.jl index ac1aa8286..96eda4c56 100644 --- a/src/odedop.jl +++ b/src/odedop.jl @@ -468,7 +468,9 @@ function odedop(coeff, p, y0, tspan; iord = 8 if h == 0.0 h = hinit(n, p, x, y, xend, posneg, k1, k2, k3, iord, hmax, abstol, reltol) - printmessages && println("hinit = $h") + if printmessages + println("hinit = $h") + end nfcn += 1 end reject = false @@ -497,10 +499,10 @@ function odedop(coeff, p, y0, tspan; F!(p, k1, x, y) end - y1, err = dopcore(coeff, n, p, x, y, h, k1, k2, k3, k4, k5, k6, k7, + y1, err, ncore = dopcore(coeff, n, p, x, y, h, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, abstol, reltol) xph = x+h - nfcn += 11 + nfcn += ncore fac11 = err^expo1 fac = fac11/facold^beta @@ -705,7 +707,7 @@ function dopcore(c::DOP853, n::Int64, p, x::Float64, y::Vector, h::Float64, k1:: deno = 1.0 end err = abs(h)*err*sqrt(1.0/(n*deno)) - return y1, maximum(abs(err)) + return y1, maximum(abs(err)), 11 end function dopcore(c::DOPRI5, n::Int64, p, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector, reltol::Vector) @@ -746,5 +748,5 @@ function dopcore(c::DOPRI5, n::Int64, p, x::Float64, y::Vector, h::Float64, k1:: err += (k4[i]/sk)*(k4[i]/sk) end err = sqrt(err/n) - return y1, maximum(abs(err)) + return y1, maximum(abs(err)), 6 end From 86832b0c77cb04ab9ef47a5f1e47a0972a38ad78 Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Fri, 27 Jun 2014 16:51:54 +0200 Subject: [PATCH 16/20] Remove magic constants. --- src/odedop.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/odedop.jl b/src/odedop.jl index 96eda4c56..f9e857306 100644 --- a/src/odedop.jl +++ b/src/odedop.jl @@ -707,7 +707,10 @@ function dopcore(c::DOP853, n::Int64, p, x::Float64, y::Vector, h::Float64, k1:: deno = 1.0 end err = abs(h)*err*sqrt(1.0/(n*deno)) - return y1, maximum(abs(err)), 11 + + # Number of function evaluations + neval = 11 + return y1, maximum(abs(err)), neval end function dopcore(c::DOPRI5, n::Int64, p, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector, reltol::Vector) @@ -748,5 +751,8 @@ function dopcore(c::DOPRI5, n::Int64, p, x::Float64, y::Vector, h::Float64, k1:: err += (k4[i]/sk)*(k4[i]/sk) end err = sqrt(err/n) - return y1, maximum(abs(err)), 6 + + # Number of function evaluations + neval = 6 + return y1, maximum(abs(err)), neval end From 134ced5286f71d8281d1bb18f5303440e30d8f82 Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Mon, 30 Jun 2014 13:47:17 +0200 Subject: [PATCH 17/20] Add stats dict, remove branch, add annotations --- src/odedop.jl | 46 +++++++++++++++++++++++----------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/src/odedop.jl b/src/odedop.jl index f9e857306..19d9fc2bb 100644 --- a/src/odedop.jl +++ b/src/odedop.jl @@ -412,12 +412,12 @@ function dopri5(f::Function, y0::Number, tspan; args...) end function odedop(coeff, p, y0, tspan; - reltol=[1e-6], abstol=[sqrt(eps())], - safe=0.9, fac1=0.333, fac2=6.0, beta=0.0, - maxstep=tspan[end]-tspan[1], initstep=0.0, - maxsteps=100000, printmessages=false, nstiff=1000, - iout=0, solout=s(x...)=return, dense=[1:length(y0)], - points=:all) + reltol::Vector{Float64}=[1e-6], abstol::Vector{Float64}=[sqrt(eps())], + safe::Float64=0.9, fac1::Float64=0.333, fac2::Float64=6.0, beta::Float64=0.0, + maxstep::Float64=tspan[end]-tspan[1], initstep=0.0, + maxsteps::Int64=100000, printmessages::Bool=false, nstiff::Int64=1000, + iout::Int64=0, solout::Function=s(x...)=return, dense::Vector{Int64}=[1:length(y0)], + points::Symbol=:all) irtrn = 0 @@ -551,13 +551,8 @@ function odedop(coeff, p, y0, tspan; for (i,j) in enumerate(dense) end end - if typeof(coeff) == DOP853 - copy!(k1, k4) - copy!(y, k5) - elseif typeof(coeff) == DOPRI5 - copy!(k1, k2) - copy!(y, y1) - end + copy!(k1, k4) + copy!(y, k5) xold = x x = xph if points == :all @@ -571,10 +566,14 @@ function odedop(coeff, p, y0, tspan; if last h = hnew idid = 1 + stats = ["function_calls"=>nfcn, + "steps"=>nstep, + "accepted"=>naccpt, + "rejected"=>nrejct] if points == :last - return x, y1 + return x, y1, stats end - return tout, yout + return tout, yout, stats end if abs(hnew) > hmax hnew = posneg*hmax @@ -727,28 +726,29 @@ function dopcore(c::DOPRI5, n::Int64, p, x::Float64, y::Vector, h::Float64, k1:: for i = 1:n y1[i] = y[i]+h*(c.a41*k1[i]+c.a42*k2[i]+c.a43*k3[i]) end - F!(p, k4, x+c.c4*h, y1) + F!(p, k8, x+c.c4*h, y1) for i = 1:n - y1[i] = y[i]+h*(c.a51*k1[i]+c.a52*k2[i]+c.a53*k3[i]+c.a54*k4[i]) + y1[i] = y[i]+h*(c.a51*k1[i]+c.a52*k2[i]+c.a53*k3[i]+c.a54*k8[i]) end - F!(p, k5, x+c.c5*h, y1) + F!(p, k9, x+c.c5*h, y1) for i = 1:n - k7[i] = y[i]+h*(c.a61*k1[i]+c.a62*k2[i]+c.a63*k3[i]+c.a64*k4[i]+c.a65*k5[i]) + k7[i] = y[i]+h*(c.a61*k1[i]+c.a62*k2[i]+c.a63*k3[i]+c.a64*k8[i]+c.a65*k9[i]) end F!(p, k6, x+h, k7) for i = 1:n - y1[i] = y[i]+h*(c.a71*k1[i]+c.a73*k3[i]+c.a74*k4[i]+c.a75*k5[i]+c.a76*k6[i]) + y1[i] = y[i]+h*(c.a71*k1[i]+c.a73*k3[i]+c.a74*k8[i]+c.a75*k9[i]+c.a76*k6[i]) end - F!(p, k2, x+h, y1) + F!(p, k4, x+h, y1) for i = 1:n - k4[i] = (c.e1*k1[i]+c.e3*k3[i]+c.e4*k4[i]+c.e5*k5[i]+c.e6*k6[i]+c.e7*k2[i])*h + k8[i] = (c.e1*k1[i]+c.e3*k3[i]+c.e4*k8[i]+c.e5*k9[i]+c.e6*k6[i]+c.e7*k4[i])*h + k5[i] = y1[i] end # Error estimation err = 0.0 for i = 1:n sk = abstol[i] + reltol[i]*max(maximum(abs(y[i])), maximum(abs(y1[i]))) - err += (k4[i]/sk)*(k4[i]/sk) + err += (k8[i]/sk)*(k8[i]/sk) end err = sqrt(err/n) From 00f4a0f5ca26b5822801d8df53b0e1b7ff2ce0f1 Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Wed, 2 Jul 2014 15:49:33 +0200 Subject: [PATCH 18/20] Fix output. --- src/odedop.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/odedop.jl b/src/odedop.jl index 19d9fc2bb..3b556a80c 100644 --- a/src/odedop.jl +++ b/src/odedop.jl @@ -465,7 +465,7 @@ function odedop(coeff, p, y0, tspan; hmax = abs(maxstep) nmax = maxsteps h = initstep - iord = 8 + iord = order(coeff) if h == 0.0 h = hinit(n, p, x, y, xend, posneg, k1, k2, k3, iord, hmax, abstol, reltol) if printmessages @@ -511,8 +511,6 @@ function odedop(coeff, p, y0, tspan; if err <= 1.0 facold = max(err, 1e-4) naccpt += 1 - F!(p, k4, xph, k5) - nfcn += 1 # Stiffness detection if mod(naccpt, nstiff) == 0 || iasti > 0 nonsti = 0 @@ -557,7 +555,7 @@ function odedop(coeff, p, y0, tspan; x = xph if points == :all push!(tout, x) - push!(yout, y1) + push!(yout, copy(k5)) end # if # solout @@ -636,6 +634,9 @@ function hinit(n::Int64, p, x::Float64, y::Vector, xend::Float64, posneg::Float6 return h*posneg end +order(c::DOP853) = 8 +order(c::DOPRI5) = 5 + function dopcore(c::DOP853, n::Int64, p, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector, reltol::Vector) y1 = 0.0 * y @@ -707,8 +708,10 @@ function dopcore(c::DOP853, n::Int64, p, x::Float64, y::Vector, h::Float64, k1:: end err = abs(h)*err*sqrt(1.0/(n*deno)) + F!(p, k4, x+h, k5) + # Number of function evaluations - neval = 11 + neval = 12 return y1, maximum(abs(err)), neval end From c0dd3b57eca97b2d8880afd65d525972204f23d2 Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Wed, 2 Jul 2014 17:57:07 +0200 Subject: [PATCH 19/20] Final numerical fixes. --- src/odedop.jl | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/src/odedop.jl b/src/odedop.jl index 3b556a80c..9611083ce 100644 --- a/src/odedop.jl +++ b/src/odedop.jl @@ -413,7 +413,8 @@ end function odedop(coeff, p, y0, tspan; reltol::Vector{Float64}=[1e-6], abstol::Vector{Float64}=[sqrt(eps())], - safe::Float64=0.9, fac1::Float64=0.333, fac2::Float64=6.0, beta::Float64=0.0, + safe::Float64=0.9, fac1::Float64=getfac1(coeff), + fac2::Float64=getfac2(coeff), beta::Float64=getbeta(coeff), maxstep::Float64=tspan[end]-tspan[1], initstep=0.0, maxsteps::Int64=100000, printmessages::Bool=false, nstiff::Int64=1000, iout::Int64=0, solout::Function=s(x...)=return, dense::Vector{Int64}=[1:length(y0)], @@ -451,7 +452,7 @@ function odedop(coeff, p, y0, tspan; copy!(y, y0) facold = 1e-4 - expo1 = 1.0/8.0 - beta*0.2 + expo1 = getexpo1(coeff, beta) facc1 = 1.0/fac1 facc2 = 1.0/fac2 posneg = sign(xend-x) @@ -569,7 +570,7 @@ function odedop(coeff, p, y0, tspan; "accepted"=>naccpt, "rejected"=>nrejct] if points == :last - return x, y1, stats + return x, k5, stats end return tout, yout, stats end @@ -636,6 +637,14 @@ end order(c::DOP853) = 8 order(c::DOPRI5) = 5 +getfac1(c::DOP853) = 0.333 +getfac1(c::DOPRI5) = 0.2 +getfac2(c::DOP853) = 6.0 +getfac2(c::DOPRI5) = 10.0 +getbeta(c::DOP853) = 0.0 +getbeta(c::DOPRI5) = 0.04 +getexpo1(c::DOP853, beta) = 1.0/8.0 - beta*0.2 +getexpo1(c::DOPRI5, beta) = 0.2 - beta*0.75 function dopcore(c::DOP853, n::Int64, p, x::Float64, y::Vector, h::Float64, k1::Vector, k2::Vector, k3::Vector, k4::Vector, k5::Vector, k6::Vector, k7::Vector, k8::Vector, k9::Vector, k10::Vector, abstol::Vector, reltol::Vector) y1 = 0.0 * y From b522eb9339b735a0a5bb3809696713cd3c5437db Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Thu, 3 Jul 2014 10:15:35 +0200 Subject: [PATCH 20/20] Remove unnecessary brackets. --- src/odedop.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/odedop.jl b/src/odedop.jl index 9611083ce..7ce15ffdf 100644 --- a/src/odedop.jl +++ b/src/odedop.jl @@ -393,7 +393,7 @@ type ODEProblemFunction f::Function end -F!(p::ODEProblemFunction, y, x, t) = copy!(y, [p.f(x, t)]) +F!(p::ODEProblemFunction, y, x, t) = copy!(y, p.f(x, t)) dop853(p, y0, tspan; args...) = odedop(dop853coeff, p, y0, tspan; args...) dopri5(p, y0, tspan; args...) = odedop(dopri5coeff, p, y0, tspan; args...)