Hi, i have to solve the nonlinear Black Scholes equation:
0=u_\tau -(sigmatilde**2/sigma**2)*u_{xx}-((sigmatilde**2/sigma**2)+D)u_x
with D=2r/sigma^2 and tau=(T/365)*(sigma^2)/2.
sigmatilde=sigma*sqrt[(1+e^{r(T-t)}a^2Ke^x*(u_xx+u_x))].
i wrote in the form:
0=u_tau -(d/dx)(Pu_x)+(d/dx)(Qu)+ Ru P=sigmatilde**2/sigma**2 #coeff of diffus.
Q=-D-sigmatilde**2/sigma**2 +(d/dx)(sigmatilde**2/sigma**2) #coeff of convection
R=(d/dx)(sigmatilde**2/sigma**2) -(d^2/dx^2)(sigmatilde**2/sigma**2) #source
The derivatives of P, Q, R are calculated numerically, Fipy representation:
- DiffusionTerm(P)+ConvectionTerm(Q)+ImplicitSourceTerm(R)+ TransientTerm()==0
The code runs but the result is not correct,
i´ve made a lot of tests and i cant
find the error, can you help me again? The code is:
L=1.
dx =0.05
nx=2.0*L/dx
mesh = Grid1D(nx = nx, dx = dx)+(-L)
r=.1,sigma=.2,T=365,K=100,S0=95,a=.02
u = CellVariable(name="Solucion numerica",mesh=mesh,value=0.)
x = mesh.getFaceCenters()[0]
Ttilde=(T/365.)*(sigma*sigma)/2.
timeStepDuration = .001
steps= Ttilde/abs(timeStepDuration)
tau=Variable()
sigmatilde = CellVariable(mesh=mesh, value=1.0)
#Initial Condition:
tmp=where( 1>numerix.exp(-x))
for i in tmp[0][:-1]:
u[i]=1.-numerix.exp(-x[i])
ux=u.getGrad()[0]
newux=ux
uxx=(u.getFaceGrad()).getDivergence()
#Here im adjusting the last value.
uxx[-1]=uxx[-2]
newuxx=uxx
if tau<2*timeStepDuration:
sigmatilde.setValue(sigma)
else:
sigmatilde.setValue(sigma*numerix.sqrt(1.+numerix.exp((2.*r*tau)/sigma**2)*a*a
*K*numerix.exp(x)[:-1]*(newux+newuxx)))
D=2*r/(sigma*sigma) ##Constant that multiplies the convective term
ratiosigmas = CellVariable(mesh=mesh, value=sigmatilde**2/sigma**2)
P = FaceVariable(mesh=mesh, value=1.0)
P=ratiosigmas #Difussion coeficient:
Q = FaceVariable(mesh=mesh, rank=1, value=-D-1.0) #Convective coeficient:
Q=(-D-ratiosigmas+ratiosigmas.getGrad())
R = CellVariable(mesh=mesh, value=1.0)
R=(ratiosigmas.getGrad())[0]-(ratiosigmas.getFaceGrad()).getDivergence()
#Boundary:
valueLeft = 0.0
valueRight = 1.-numerix.exp(-D*tau-x)
BCs = (FixedValue(faces=mesh.getFacesRight(), value=valueRight),
FixedValue(faces=mesh.getFacesLeft(), value=valueLeft))
dt=timeStepDuration
while tau()< Ttilde:
tau.setValue(tau()+dt)
uxx=(u.getFaceGrad()).getDivergence()
uxx[-1]=uxx[-2]
newuxx=uxx
if tau<2.*timeStepDuration:
sigmatilde.setValue(sigma)
else:
sigmatilde.setValue(sigma*numerix.sqrt(1.+numerix.exp((2.*r*tau)/sigma**2)
*a**2*K*numerix.exp(x)[:-1]*(newux+newuxx)))
ratiosigmas=sigmatilde**2/sigma**2
P=ratiosigmas
Q=(-D-ratiosigmas+ratiosigmas.getGrad())
R=(ratiosigmas.getGrad())[0]-(ratiosigmas.getFaceGrad()).getDivergence()
eqX.solve(var=u,boundaryConditions=BCs,dt=timeStepDuration)
if __name__ == '__main__':
viewer.plot()
