In [1]:
using Plots,LaTeXStrings
default(markersize=3,linewidth=1.5)
using LinearAlgebra,DifferentialEquations
using SparseArrays
include("FNC.jl");
Here are the parameters for a particular simulation of the Black--Scholes equation.
In [2]:
Smax = 8; T = 6;
K = 3; sigma = 0.06; r = 0.08;
We discretize space and time.
In [3]:
m = 200; h = Smax / m;
x = h*(0:m)
n = 1000; tau = T / n;
t = tau*(0:n)
lambda = tau / h^2; mu = tau / h;
Set the initial condition, then march forward in time.
In [4]:
V = zeros(m+1,n+1)
V[:,1] = @. max( 0, x-K )
for j = 1:n
# Fictitious value from Neumann condition.
Vfict = 2*h + V[m,j]
Vj = [ V[:,j]; Vfict ]
# First row is zero by the Dirichlet condition.
for i = 2:m+1
diff1 = (Vj[i+1] - Vj[i-1])
diff2 = (Vj[i+1] - 2*Vj[i] + Vj[i-1])
V[i,j+1] = Vj[i] +
(lambda*sigma^2*x[i]^2/2)*diff2 +
(r*x[i]*mu)/2*diff1 - r*tau*Vj[i]
end
end
We plot at a few times.
In [5]:
select_times = @. 1 + 250*(0:4)
show_times = t[select_times]
plot(x,V[:,select_times],label=["t=$t" for t in show_times],
title="Black-Scholes solution",leg=:topleft,
xaxis=("stock price"),yaxis=("option value") )
Out[5]:
The lowest curve is the initial condition, and the highest curve is the last time. The results are easy to interpret, recalling that the time variable really means "time before strike." Say you are close to the option's strike time. If the current stock price is, say, $S=2$, then it's not likely that the stock will end up over the strike price $K=3$ and therefore the option has little value. On the other hand, if presently $S=3$, then there are good odds that the option will be exercised at the strike time, and you will need to pay a substantial portion of the stock price in order to take advantage.
Let's try to extend the simulation time to $T=8$, keeping everything else the same.
In [6]:
T = 8;
n = 1000; tau = T / n;
t = tau*(0:n)
lambda = tau / h^2; mu = tau / h;
for j = 1:n
# Fictitious value from Neumann condition.
Vfict = 2*h + V[m,j]
Vj = [ V[:,j]; Vfict ]
# First row is zero by the Dirichlet condition.
for i = 2:m+1
diff1 = (Vj[i+1] - Vj[i-1])
diff2 = (Vj[i+1] - 2*Vj[i] + Vj[i-1])
V[i,j+1] = Vj[i] +
(lambda*sigma^2*x[i]^2/2)*diff2 +
(r*x[i]*mu)/2*diff1 - r*tau*Vj[i]
end
end
plot(x,V[:,select_times],label=["t=$t" for t in show_times],
title="Black-Scholes solution",leg=:topleft,
xaxis=("stock price"),yaxis=("option value",[0,1e20]) )
Out[6]:
This "solution" is nonsensical. Look at the scale of the ordinate!
Let's try out the Euler and backward Euler time stepping methods using the second-order semidiscretization:
In [7]:
m = 100
x,Dx,Dxx = FNC.diffper(m,[0,1]);
First we apply the Euler discretization.
In [8]:
tfinal = 0.05; n = 500;
tau = tfinal/n; t = tau*(0:n);
U = zeros(m,n+1);
This is where we set the initial condition. It isn't mathematically periodic, but the end values and derivatives are so small that for numerical purposes it may as well be.
In [9]:
U[:,1] = @. exp( -60*(x-0.5)^2 );
The Euler time stepping simply multiplies by a constant matrix for each time step.
In [10]:
A = I + tau*Dxx;
for j = 1:n
U[:,j+1] = A*U[:,j]
end
Things seem to start well.
In [11]:
plot(x,U[:,1:3:7],label=["t=$t" for t in t[1:3:7]],
xaxis=(L"x"),yaxis=(L"u(x,t)"),
title="Heat equation by forward Euler")
Out[11]:
Shortly thereafter, though, there is nonphysical growth.
In [12]:
plot(x,U[:,13:3:19],label=["t=$t" for t in t[13:3:19]],
xaxis=(L"x"),yaxis=(L"u(x,t)"),
title="Heat equation by forward Euler")
Out[12]:
The growth is exponential in time.
In [13]:
M = vec(maximum( abs.(U), dims=1 )) # max in each column
plot(t,M,xaxis=(L"t"), yaxis=(:log10,L"\max_x |u(x,t)|"),
title="Nonphysical growth",leg=:none)
Out[13]:
Now we try backward Euler. In this case there is a tridiagonal linear system to solve at each time step. We will use a sparse matrix to get sparse LU factorization, although the time savings at this size are negligible.
In [14]:
B = sparse(I - tau*Dxx)
for j = 1:n
U[:,j+1] = B\U[:,j]
end
plot(x,U[:,1:100:501],label=["t=$t" for t in t[1:100:501]],
xaxis=(L"x"),yaxis=(L"u(x,t)"),
title="Heat equation by backward Euler")
Out[14]:
This solution looks physically realistic, as the large concentration in the center diffuses outward. Observe that the solution remains periodic in space.
We set up the semidiscretization and initial condition in $x$ just as before.
In [15]:
m = 100
x,Dx,Dxx = FNC.diffper(m,[0,1])
u0 = @. exp( -60*(x-0.5)^2 );
Now, though, we apply a DifferentialEquations solver to the initial-value problem $\mathbf{u}'=\mathbf{D}_{xx}\mathbf{u}$.
In [16]:
tfinal = 0.05
ODE = (u,p,t) -> Dxx*u;
IVP = ODEProblem(ODE,u0,(0,tfinal))
u = solve(IVP,DP5()) # select DP5 solver
plot(x,u[:,1:50:end],label=["t=$t" for t in round.(u.t[1:50:end],digits=4)],
xaxis=(L"x"), yaxis=(L"u(x,t)"),
title="Heat equation by ode45")
Out[16]:
The solution is at least plausible. But the number of time steps that were selected automatically is surprisingly large, considering how smoothly the solution changes.
In [17]:
num_steps_1 = length(u.t)-1
Out[17]:
Now we apply another solver.
In [18]:
u = solve(IVP,Rodas4P());
The number of steps selected is reduced by a factor of more than 40!
In [19]:
num_steps_2 = length(u.t)-1
Out[19]:
Both time stepping methods solved $\mathbf{u}'=\mathbf{D}_{xx}\mathbf{u}$, for the matrix
In [20]:
m = 40; Dxx = FNC.diffper(m,[0,1])[3];
The eigenvalues of this martrix are real and negative.
In [21]:
lambda = eigvals(Dxx)
scatter(real(lambda),imag(lambda),
xaxis=("Re \\lambda"), yaxis=("Im \\lambda"),
title="Eigenvalues",leg=:none)
Out[21]:
The Euler method is absolutely stable in the region $|\zeta+1| \le 1$ in the complex plane:
In [22]:
phi = 2pi*(0:360)/360
z = @.exp(1im*phi) - 1; # unit circle shifted to the left by 1
plot(Shape(real(z),imag(z)),color=RGB(.8,.8,1),
xaxis=("Re \\lambda"), yaxis=("Im \\lambda"), aspect_ratio=1,
title="Stability region",leg=:none)
Out[22]:
In order to get inside this region, we have to find $\tau$ such that $\lambda \tau > -2$ for all eigenvalues $\lambda$. This is an upper bound on $\tau$.
In [23]:
lambda_min = minimum(lambda)
@show max_tau = -2 / lambda_min;
Here we plot the resulting values of $\zeta=\lambda \tau$.
In [24]:
zeta = lambda*max_tau
scatter!(real(zeta),imag(zeta),
title="Stability region and \\zeta values")
Out[24]:
In backward Euler, the region is $|\zeta-1|\ge 1$. Because they are all on the negative real axis, all of the $\zeta$ values will fit no matter what $\tau$ is chosen.
In [25]:
plot(Shape([-6,6,6,-6],[-6,-6,6,6]),color=RGB(.8,.8,1))
z = @.exp(1im*phi) + 1; # unit circle shifted right by 1
plot!(Shape(real(z),imag(z)),color=:white)
scatter!(real(zeta),imag(zeta),
xaxis=([-4,2],"Re \\lambda"), yaxis=([-3,3],"Im \\lambda"),aspect_ratio=1,
title="Stability region and \\zeta values",leg=:none)
Out[25]:
In [26]:
f = function (u,p,t)
q,s,w = p
[ s*(u[2]-u[1]*u[2]+u[1]-q*u[1]^2);
(-u[2]-u[1]*u[2]+u[3])/s;
w*(u[1]-u[3]) ]
end
Out[26]:
One of the stiff solvers, called Rodas4P, is fast.
In [27]:
IVP = ODEProblem(f,[1.,1.,4.],(0.,6.),(8.375e-6,77.27,0.161))
@elapsed u = solve(IVP,Rodas4P()) # select Rodas4P solver
Out[27]:
In [28]:
plot(u)
scatter!(u.t,u[:,:]',m=2,xaxis=(L"t"),yaxis=(L"u(t)"),
title="Oregonator",leg=:none)
Out[28]:
In [29]:
@show num_steps_stiff = length(u.t);
You can see that the stiff solver takes small time steps only when some part of the solution is changing rapidly. However, a nonstiff solver, based on an explicit method, is much slower and takes a lot more steps.
In [30]:
@elapsed u = solve(IVP,DP5()) # select DP5 solver
Out[30]:
In [31]:
@show num_steps_nonstiff = length(u.t);
Here is the Jacobian matrix at any value of $\mathbf{u}$.
In [32]:
q,s,w = (8.375e-6,77.27,0.161)
J = u -> [ -s*(u[2]+1-2*q*u[1]) s*(1-u[1]) 0;
-u[2]/s (-1-u[1])/s 1/s;
w 0 -w];
During the early phase, the time steps seem fairly large. The eigenvalues around $t=1/2$ are
In [33]:
t = u.t
i1 = findfirst(@. t>0.5)
lambda1 = eigvals( J(u[:,i1]) )
Out[33]:
These are real and negative. Checking the stability region of RK4 along the negative real axis, we see that stability requires a maximum time step
In [34]:
maxstep1 = 2.8 / maximum(abs.(lambda1))
Out[34]:
The actual step size chosen by the solver was comparable:
In [35]:
step1 = t[i1+1] - t[i1]
Out[35]:
Later in the simulation, the steps seem quite small compared to the apparent rate of activity. We look near $t=4$:
In [36]:
i2 = findfirst(@. t>4)
lambda2 = eigvals( J(u[:,i2]) )
Out[36]:
These are also real and negative. We compare the maximum and observed step sizes again:
In [37]:
@show maxstep2 = 2.8 / maximum(abs.(lambda2));
@show step2 = t[i2+1] - t[i2];
Again, the eigenvalues give a good indication of how the steps are being chosen, at least to order of magnitude.
We solve $u_t=u_{xx}$ on $[-1,1]$ subject to the Dirichlet conditions $u(-1,t)=0$, $u(1,t)=2$.
In [38]:
m = 100
x,Dx,Dxx = FNC.diffcheb(m,[-1,1]);
Our next step is to write a function that defines $\mathbf{f}$. Since the boundary values are given explicitly, there is no need to "solve" for them---we just append them to each end of the vector.
In [39]:
extend(v) = [0;v;2];
We can also define the inverse operation of chopping off the boundary values from a full vector. Note that the indexing starts at 1 and ends at $m+1$ for the extended vector.
In [40]:
chop(u) = u[2:m];
All the pieces are now in place to define and solve the IVP.
In [41]:
f = function (v,p,t)
u = extend(v)
phi = Dxx*u
return chop(phi)
end
u0 = @. 1 + sin(pi/2*x) + 3*(1-x^2)*exp(-4*x^2);
V = solve(ODEProblem(f,chop(u0),(0.,0.15)));
We extend the solution to the boundaries at each time, then plot.
In [42]:
t = range(0,stop=0.15,length=11)
U = hcat( (extend(V(t)) for t in t)... )
plot(x,U,label=["t=$t" for t in t],
xaxis=(L"x"), yaxis=(L"u(x,t)"), title="Heat equation",leg=:topleft)
Out[42]:
We solve a diffusion equation with source term: $u_t=u^2+u_{xx}$, on $[-1,1]$ subject to homogeneous Dirichlet conditions.
In [43]:
m = 100
x,Dx,Dxx = FNC.diffcheb(m,[-1,1]);
In [44]:
extend(v) = [0;v;0]; # extend to boundary
chop(u) = u[2:m]; # discard boundary
ODE = function (v,p,t)
u = extend(v)
uxx = Dxx*u
f = @.u^2 + uxx
return chop(f)
end
Out[44]:
All the pieces are now in place to define and solve the IVP.
In [45]:
u0 = @. 6*(1-x^2)*exp(-4*(x-.5)^2)
V = solve(ODEProblem(ODE,chop(u0),(0.,1.5)));
Extend the solution to the boundaries at each time, then plot.
In [46]:
t = range(0,stop=1.5,length=7)
U = hcat( (extend(V(t)) for t in t)... )
plot(x,U,label=["t=$t" for t in t],
xaxis=(L"x"), yaxis=(L"u(x,t)"), title="Heat equation with source",leg=:topleft)
Out[46]:
An animation better captures how the source term takes over gradually but at an accelerating rate.
In [47]:
an = @animate for t = range(0,stop=1.5,length=100)
plot(x,extend(V(t)),label="t=$(round(t,digits=3))",
xaxis=(L"x"), yaxis=([0,10],L"u(x,t)"), title="Heat equation with source",leg=:topleft)
end
gif(an,"molbratu.gif")
Out[47]:
In [48]:
Smax = 8; T = 8;
K = 3; sigma = 0.06; r = 0.08;
m = 200; x,Dx,Dxx = FNC.diffmat2(m,[0,Smax]);
h = x[2]-x[1];
Using the boundary conditions and defining the ODE follow next.
In [49]:
extend(v) = [ 0; v; 2/3*(h-0.5*v[m-2]+2*v[m-1]) ];
chop(u) = u[2:m];
ODE = function (v,p,t)
sigma,r = p
u = extend(v)
ux = Dx*u; uxx = Dxx*u;
f = @.sigma^2/2*x^2*uxx + r*x*ux - r*u
return chop(f)
end;
Now we define the initial conditions and solve the IVP.
In [50]:
v0 = @.max( 0, x-K )
V = solve(ODEProblem(ODE,chop(v0),(0.,T),(sigma,r)));
Extend the solution to the boundaries at each time, then plot.
In [51]:
t = range(0,stop=T,length=7)
U = hcat( (extend(V(t)) for t in t)... )
plot(x,U,label=["t=$(round(t,digits=2))" for t in t],
xaxis=(L"x"), yaxis=(L"u(x,t)"), title="Black-Scholes equation",leg=:topleft)
Out[51]:
In [52]:
an = @animate for t = range(0,stop=T,length=100)
plot(x,extend(V(t)),label="t=$(round(t,digits=3))",
xaxis=(L"x"), yaxis=([0,10],L"u(x,t)"), title="Black-Scholes equation",leg=:topleft)
end
gif(an,"molbsbc.gif")
Out[52]: