source('./improved_euler.R')

t0 = 0
y0 = 1
t1 = 0.4
h = 0.1

# Where do we want to evaluate the solution:
#
x_grid = seq( t0, t1, by=h )

# Problem 16:
#
# The exact solution:
#
y_exact_fn = function(t){ (1 + exp(2*t))/2 }
y_exact = y_exact_fn(x_grid)

# Compute the local truncation error bound for the improved Euler method: (2/3) e^{2 \bar{t}_n} h^3: 
#
lteb = (2/3) * exp( 2*(x_grid + h) ) * h^3 # holds |e_1), |e_2|, |e_3|, |e_4|,
lteb = lteb[c(1, 4)]

# Compute the two true errors at y_1 and y_4:
#
dy.dt = function(t, y){ 2*y - 1 }

res = improved_euler(dy.dt, h=h, start=t0, y0=y0, end=t0+h)
e_y1 = res$ys[2]

t3 = t0 + 3*h
y3 = y_exact_fn(t3)
res = improved_euler(dy.dt, h=h, start=t3, y0=y3, end=t3+h)
e_y4 = res$ys[2]

# The true error for that improved Euler step:
#
true_error = y_exact[c(1, 4)+1] - c( e_y1, e_y4 )

A = rbind( lteb, abs(true_error) )
rownames(A) = c('Local Truncation Error:', 'abs(True Error):')
colnames(A) = c( 'abs(e1)', 'abs(e4)' )
print(A)


# Problem 17:
#

# The exact solution:
#
y_exact_fn = function(t){ t/2 + exp(2*t) }
y_exact = y_exact_fn(x_grid)

# Compute the local truncation error bound for the improved Euler method: (4/3) e^{2 \bar{t}_n} h^3
#
lteb = (4/3) * exp( 2*(x_grid + h) ) * h^3 # holds le_l|, |e_2|, |e_3|, |e_4|, ...
lteb = lteb[c(1, 4)]

# Compute the two true errors at y_1 and y_4:
#
dy.dt = function(t, y){ 1/2 - t + 2*y }

res = improved_euler(dy.dt, h=h, start=t0, y0=y0, end=t0+h)
e_y1 = res$ys[2]

t3 = t0 + 3*h
y3 = y_exact_fn(t3)

res = improved_euler(dy.dt, h=h, start=t3, y0=y3, end=t3+h)
e_y4 = res$ys[2]

# The true error for that improved Euler step:
#
true_error = y_exact[c(1, 4)+1] - c( e_y1, e_y4 )

A = rbind( lteb, abs(true_error) )
rownames(A) = c('Local Truncation Error:', 'abs(True Error):')
colnames(A) = c( 'abs(e1)', 'abs(e4)' )
print(A)