1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
|
// file BlackScholes2D.edp
int m=30,L=80,LL=80, j=100;
real sigmax=0.3, sigmay=0.3, rho=0.3, r=0.05, K=40, dt=0.01;
mesh th=square(m,m,[L*x,LL*y]);
fespace Vh(th,P1);
Vh u=max(K-max(x,y),0.);
Vh xveloc, yveloc, v,uold;
for (int n=0; n*dt <= 1.0; n++)
{
if(j>20) { th = adaptmesh(th,u,verbosity=1,abserror=1,nbjacoby=2,
err=0.001, nbvx=5000, omega=1.8, ratio=1.8, nbsmooth=3,
splitpbedge=1, maxsubdiv=5,rescaling=1) ;
j=0;
xveloc = -x*r+x*sigmax^2+x*rho*sigmax*sigmay/2;
yveloc = -y*r+y*sigmay^2+y*rho*sigmax*sigmay/2;
u=u;
};
uold=u;
solve eq1(u,v,init=j,solver=LU) = int2d(th)( u*v*(r+1/dt)
+ dx(u)*dx(v)*(x*sigmax)^2/2 + dy(u)*dy(v)*(y*sigmay)^2/2
+ dy(u)*dx(v)*rho*sigmax*sigmay*x*y/2 + dx(u)*dy(v)*rho*sigmax*sigmay*x*y/2)
+ int2d(th)( -v*convect([xveloc,yveloc],dt,uold)/dt)+ on(2,3,u=0);
j=j+1;
};
plot(u,wait=1,value=1);
plot(th,wait=1);
|
