# Stochastic Calculus, Courant Institute, NYU, Fall 2013 # http://www.math.nyu.edu/faculty/goodman/teaching/StochCalc2013/ # Assignment 6, Solve the backward equation and compare to Monte Carlo evaluation # (your name and email here) #---------------------------------------------------------------------------- # The possible value functions tooth <- function( x, a){ # Evaluate a function of x defined for -a <= x <= a. # In this case, the sawtooth function, V(x) = 0 if x < 0, V(x) = 1-x/a if x > 0 if ( x > 0 ) return (1.-(x/a) ) return( 0.) } cosf <- function( x, a){ # Evaluate a function of x defined for -a <= x <= a. # In this case, the cosine function, V(x) = cos(pi*x/(2*a) return( cos(3.1415926535787*x/(2*a) ) ) } cosSol <- function( xp, a, T, N){ # The exact solution for the cosine final condition fe = array( 0, N) k = 3.1415926535787/(2*a) for ( j in 1:N ){ fe[j] = cos( k*xp[j] ) * exp( -k*k*T/2) } return( fe ) } #---------------------------------------------------------------------------- # A procedure that sets up a grid, mostly for plotting return an array of the x values xvals <- function(a, N){ # See below for the definition of a and N dx = 2*a/(N-1) x = -a # x moves from left to right, starting at -a xp = array( 0, N) for ( j in 1:N ){ xp[j] = x # Store the x values, for plotting x = x + dx # move to the next grid point } return(xp) } #----------------------------------------------------------------------------- # The code that solves the backward equation using finite differences fds <- function( T, # The final time a, # The half-width of the interval N, # The number of grid points in space, including the boundary points xp, # The x values of the grid points V){ # A function that evaluates the final condition: f(x,T) = V(x,a), a = interval half-width dx = 2*a/(N-1) # The distance between grid points dt = .5*dx*dx # The time step fnt = T/dt # A floating point number giving the number of time steps, probably not an integer nt = as.integer(fnt) + 1 # Round up to the nearest integer. The as.integer function rounds down t = T/nt # Adjust the time step (a little) to be consistent with nt ft = array( 0 ,N) # The finite difference value function at time t <= T ftm1 = array( 0, N) # The finite difference value function at time t - dt # Set the final condition for ( j in (1:N)) { ft[j] = V(xp[j],a) # The final condition, as a function of x } ft[1] = 0. # boundary conditions at the left ... ft[N] = 0. # ... and the right end. # Compute the finite difference coefficients f_{t-dt,j} = cp*f_{t,j+1} + c0*f_{t,j} + cm*f_{t,j-1} cp = dt/(2*dx*dx) c0 = 1. - (dt/(dx*dx)) cm = dt/(2*dx*dx) for ( k in (1:nt) ) { # nt is the number of time steps for ( j in (2:(N-1) )) { # loop over the non-boundary points in x ftm1[j] = cp*ft[j+1] + c0*ft[j] + cm*ft[j-1] } ftm1[1] = 0. # put in the boundary values ... ftm1[N] = 0. # ... do not have to do this every time. It is here for clearity. for (j in (1:N) ){ # copy ftm1 back to ft as t is replaced by t - dt ft[j] = ftm1[j] } } return(ft) } #----------------------------------------------------------------------------- # The main program # Set parameters T = 10. # The final time (needs a decimal point) a = 3. # An absorbing boundary at -a and a N = 12 # The number grid points in [-a, a] #---------------------------------------------------------------------------- # Experiment 1: Cosine final condition, compare numerical and exact solution xp = xvals( a, N) # All the x values, for plotting ft = fds( T, a, N, xp, cosf) # Calculate the finite difference solution on a coarse grid Nf = 2*N - 1 # The number of points for the fine grid, to make twice as many intervals xpf = xvals( a, Nf) ftf = fds( T, a, Nf, xpf, cosf) # Calculate the finite difference solution on a finer grid fe = cosSol( xp, a, T, N) setwd("/Users/jg/desktop/notes/StochCalc2013/codes") # set the directory pdf( "Assignment6Figure1.pdf") # open and name the plot file title = sprintf('Cosine final condition, difference and exact solution') PlotInfo = sprintf('T = %8.2f, N = %6d', T, N) xlabel = "x" ylabel = "f" ApproxCh = 20 # Character for approximate solution points ExactCh = 10 # Character for exact solution points plot( xp, ft, # plot the computed solution main = title, # overall plot title sub = PlotInfo, # subtitle has more info on the plot xlab = xlabel, # label for the x axis ylab = ylabel, # label for the y axis type = 'p', # draw symbols pch = ApproxCh, # a kind of box col = 'blue' ) # blue points (xp, fe, pch = ExactCh) LegendText = c("approx soln","exact soln") legend( "topleft", LegendText, pch = c(ApproxCh,ExactCh)) dev.off( which = dev.cur() ) # This closes the current "device", which ... Nf = 2*N - 1 # The number of points for the fine grid, to make twice as many intervals xpf = xvals( a, Nf) ftf = fds( T, a, Nf, xpf, cosf) # Calculate the finite difference solution on a finer grid fe = cosSol( xp, a, T, N) #---------------------------------------------------------------------------- # Experiment 2: Hat function final condition, compare coarse and fine grid numerics T = .1 N = 50 Nf = 2*N - 1 # The number of points for the fine grid, to make twice as many intervals xp = xvals( a, N) # All the x values, for plotting xpf = xvals( a, Nf) # The fine grid points ft = fds( T, a, N , xp , tooth) # Calculate the finite difference solution on a original grid ftf = fds( T, a, Nf, xpf, tooth) # Calculate the finite difference solution on a finer grid pdf( "Assignment6Figure2.pdf") # open and name the plot file title = sprintf('Sawtooth final condition, coarse and fine approx') PlotInfo = sprintf('T = %8.2f, N = %6d and N = %6d', T, N, Nf) xlabel = "x" ylabel = "f" CoarseCh = 20 # Character for the coarse grid approximation FineCh = 10 # Character for the fine grid approximation plot( xpf, ftf, # plot the computed solution main = title, # overall plot title sub = PlotInfo, # subtitle has more info on the plot xlab = xlabel, # label for the x axis ylab = ylabel, # label for the y axis type = 'p', # draw symbols pch = CoarseCh, # a kind of box col = 'blue' ) # blue points (xp, ft, pch = FineCh) LegendText = c("fine grid soln","coarse grid soln") legend( "topleft", LegendText, pch = c(CoarseCh,FineCh), col = c('blue', 'black')) dev.off( which = dev.cur() ) # This closes the current "device", which ... #---------------------------------------------------------------------------- # Experiment 3: Hat function final condition, compare solutions at different times T1 = .02 # Final times for the various solutions T2 = .2 T3 = 2. T4 = 10. N = 100 xp = xvals( a, N) # All the x values, for plotting ft1 = fds( T1, a, N, xp, tooth) # Calculate the finite difference solutions ft2 = fds( T2, a, N, xp, tooth) ft3 = fds( T3, a, N, xp, tooth) ft4 = fds( T4, a, N, xp, tooth) pdf( "Assignment6Figure3.pdf") # open and name the plot file title = sprintf('Sawtooth final condition, several solutions') PlotInfo = sprintf('T1 =%8.2f, T2 =%8.2f, T3 =%8.2f, T4 =%8.2f', T1, T2, T3, T4) Col1 = 'blue' # The colors of the various solutions Col2 = 'red' Col3 = 'green' Col4 = 'black' xlabel = "x" ylabel = "f" PlotCh = 20 # Character for approximate solution points, a dot plot( xp, ft1, # plot the T1 solution main = title, # overall plot title sub = PlotInfo, # subtitle has more info on the plot xlab = xlabel, # label for the x axis ylab = ylabel, # label for the y axis type = 'p', # draw symbols pch = PlotCh, # plot symbol col = Col1 ) # solution colors points (xp, ft2, pch = PlotCh, col = Col2) # Add the other solutions points (xp, ft3, pch = PlotCh, col = Col3) points (xp, ft4, pch = PlotCh, col = Col4) LegendText = c("Time 1 soln","Time 2 soln","Time 3 soln","Time 4 soln") legend( "topleft", # Where the legend goes LegendText, # duh pch = c(PlotCh,PlotCh,PlotCh,PlotCh), # the plot characters, all the same here col = c(Col1, Col2, Col3, Col4)) # the solution colors dev.off( which = dev.cur() ) # This closes the current "device", which ...