Skip to content

Commit

Permalink
a working, though incorrect, Netwon iteration
Browse files Browse the repository at this point in the history
  • Loading branch information
@obiajulu
obiajulu committed Jul 8, 2016
1 parent 0a39b32 commit b8a898c
Showing 1 changed file with 53 additions and 42 deletions.
95 changes: 53 additions & 42 deletions src/radau.jl
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
#=
(1) done()
✓(1) done()
(2) trialstep!()
(3) errorcontol!()
(4) odercontrol!()
Expand All @@ -14,6 +14,7 @@ using Polynomials
using ForwardDiff
using ODE
using Parameters
using Compat


immutable TableauRKImplicit{Name, S, T} <: ODE.Tableau{Name, S, T}
Expand Down Expand Up @@ -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 ],
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -108,46 +114,43 @@ 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
else
rollback!()
end
return status()
return status()=#
end
end

Expand All @@ -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
Expand All @@ -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
Expand All @@ -204,35 +212,26 @@ 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

# 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)
Expand All @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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

0 comments on commit b8a898c

Please sign in to comment.