diff --git a/src/radau.jl b/src/radau.jl index e2bd2c2b4..a9d282441 100644 --- a/src/radau.jl +++ b/src/radau.jl @@ -1,5 +1,5 @@ #= - (1) done() + ✓(1) done() (2) trialstep!() (3) errorcontol!() (4) odercontrol!() @@ -14,6 +14,7 @@ using Polynomials using ForwardDiff using ODE using Parameters +using Compat immutable TableauRKImplicit{Name, S, T} <: ODE.Tableau{Name, S, T} @@ -44,13 +45,13 @@ function TableauRKImplicit(name::Symbol, order::Integer, T::Type, end ## Tableaus for implicit RK methods -const bt_radau3 = TableauRKImplicit(:radau3,3, Rational{Int64}, +const bt_radau3 = TableauRKImplicit(:radau3,3, Float64, [5//12 -1//12 3//4 1//4], [3//4, 1//4]', [1//3, 1]) -const bt_radau5 = TableauRKImplicit(:radau5,5, Rational{Int64}, +const bt_radau5 = TableauRKImplicit(:radau5,5, Float64, [11/45-7*√(6)/360 37/225-169*√(6)/1800 -2/225 + √(6)/75 37/225+169*√(6)/1800 11/45+7*√(6)/360 -2/225 - √(6)/75 4/9-√(6)/36 4/9+√(6)/36 1/9 ], @@ -62,6 +63,7 @@ const bt_radau5 = TableauRKImplicit(:radau5,5, Rational{Int64}, # State for Radau Solver ########################################### type RadauState{T,Y} + f:: Function tfinal ::T tdir :: Integer minstep ::Float64 @@ -82,6 +84,8 @@ type RadauState{T,Y} btab #<:TableauRKImplicit # tableau according to stage number stageNum ::Integer M ::Array{Float64,2} + reltol::Float64 + abstol::Float64 # work arrays end @@ -90,7 +94,9 @@ end # Radau Solver ########################################### function ode_radau(f, y0, tspan, order ::Integer = 5; M = eye(length(y0),length(y0)), - minstep = abs(tspan[length(tspan)]-tspan[1])/10^12 + minstep = abs(tspan[length(tspan)]-tspan[1])/10^12, + reltol = 1.0e-5, + abstol = 1.0e-8, ) # Set up stageNum = Integer((order+1)/2) @@ -108,38 +114,35 @@ function ode_radau(f, y0, tspan, order ::Integer = 5; M = eye(length(y0),lengt dy = f(t, y) ## previous data set to null to begin - tpre = NaN - ypre = zeros(y0) - dypre = zeros(y0) + tpre = tspan[1] + ypre = deepcopy(y0) + dypre = f(t, y) step = 1 finished = false - stageNum = stageNum # get right tableau for stage number - if stageNum ==3 + if order ==3 btab = bt_radau3 - elseif stageNum ==5 + elseif order ==5 btab = bt_radau5 else - @show typeof(stageNum) - @show constRadauTableau(stageNum) btab = constRadauTableau(stageNum) end - @show typeof(tfinal) ## intialize state - st = RadauState{T,Y}(tfinal,tdir, minstep,h, + st = RadauState{T,Y}(f, tfinal,tdir, minstep,h, t, y, dy, tpre, ypre, dypre, step, finished, - btab, stageNum, M) + btab, stageNum, M, + reltol, abstol) - println("Good so far!") # Time stepping loop while !done_radau(st) + @show st.step stats = trialstep!(st) - err, stats, st = errorcontrol!(st) # (1) adaptive stepsize (2) error + #=err, stats, st = errorcontrol!(st) # (1) adaptive stepsize (2) error if err < 1 stats, st = ordercontrol!() accept = true @@ -147,7 +150,7 @@ function ode_radau(f, y0, tspan, order ::Integer = 5; M = eye(length(y0),lengt rollback!() end - return status() + return status()=# end end @@ -171,11 +174,11 @@ end # Time vector going from 0 to 2*PI in 0.01 increments t = 0:0.1:4*pi; -ode_radau(f,y,t,7) +ode_radau(f,y,t,5) ########################################### # iterator like helper functions ########################################### - +"done function for when to stop radau time stepping loop" function done_radau(st) @unpack st: h, t, tfinal, minstep, tdir if h < minstep || tdir*t >= tdir*tfinal @@ -185,14 +188,19 @@ function done_radau(st) end end -"trial step for ODE with mass matrix" +"trial step for ODE with mass matrix with Radau IIA solver" function trialstep!(st) - @unpack st: h, t,y, tfinal, btab - + @unpack st: h, t,y,dy, tfinal, btab, f, stageNum,M,reltol,abstol, ypre + dof = length(y) # Form ⃗z from y - - # Calculate Jacobian - g(z) = f(t,z) + z = zeros(stageNum) #TODO: use better initial values for z + zpre = zeros(stageNum) + Δzpre = zeros(stageNum) + Δz = zeros(stageNum) + κ = .01 # See Harier Vol II pg 121 + + # Calculate Jacobian at (t_n, y_n) + g(y_val) = f(t,y_val) J = ForwardDiff.jacobian(g, y) # Use to form simplied Netwon iteration matrix @@ -204,7 +212,6 @@ function trialstep!(st) # I_N = eye(dof,dof) I_s = eye(stageNum,stageNum) - M = rand(dof,dof) AoplusJ = kron(btab.a,J) IoplusM = kron(I_s,M) G = IoplusM-h*AoplusJ @@ -212,27 +219,19 @@ function trialstep!(st) # Use Netwon interation (TODO: use the transformation T^(-1)A^(-1)T = Λ, # W^k = (T^-1 ⊕ I)Z^k version of iteration) - ## initial variables iteration - ## Initialize - z = zeros(dof) #TODO: use better initial values for z - zpre = zeros(dof) - Δzpre = zeros(dof) - Δz = zeros(dof) - κ - ## Matrices used for one round of iteration Ginv = inv(G) - Ginv_block = Array{Float64,2}[Ginv[i*stageNum + [1:dof], j*stageNum+[1:dof]] for i = 0:stageNum-1, j= 0:stageNum-1] - AoplusI_block = Array{Float64,2}[btab.a[i,j]*I_N for i=1:stageNum, j=1:stageNum] + @compat Ginv_block = Array{Float64,2}[Ginv[i*dof + [1:dof], j*dof+[1:dof]] for i = 0:stageNum-1, j= 0:stageNum-1] + @compat AoplusI_block = Array{Float64,2}[btab.a[i,j]*I_N for i=1:stageNum, j=1:stageNum] iterate = true count = 0 - while iterate - Δz = Ginv_block*(-zpre + h*AoplusI_block*F(f,z,y,t,c,h)) + while iterate && count<=7 + Δz = Ginv_block*(-zpre + h*AoplusI_block*F(f,z,y,t,btab.c,h)) z = zpre + Δz # ⃗z^(k+1) = ⃗z ^ (k) - Δ⃗z^(k) # Stop condition for the Netwon iteration - if (count >=1) + if (count >=2) Θ = norm(Δz)/norm(Δzpre) η = Θ/(1-Θ) if η*norm(Δz) <= κ*min(reltol,abstol) @@ -244,7 +243,7 @@ function trialstep!(st) zpre = z Δzpre = Δz - count+=1 + @show count+=1 end # Once Netwon method converges after some m steps, estimated next step size @@ -256,6 +255,18 @@ function trialstep!(st) for i = 1 : stageNum ynext += z[i]*d[i] end + + #update state + st.ypre = y + st.dypre = dy + st.tpre = t + + st.y = ynext + st.t = t+h + st.dy = f(t+h,ynext) + + st.step = st.step+1 + return (count) end function errorcontrol!(st) @@ -363,5 +374,5 @@ end " Calculates the array of derivative values between t and tnext" function F(f,z,y,t,c,h) - return Array{Float64,1}[f(t+c[i]*h, y+z[i]) for i=1:stageNum] + return Array{Float64,1}[f(t+c[i]*h, y+z[i]) for i=1:length(z)] end