Dear FreeFem++ Community,

I am a Finance student exploring the partial differential equation (PDE) approach for pricing American basket options using FreeFem++. I came across a script named “american.edp” and am seeking clarification on a specific code segment that appears to deviate from the theory referenced in a presentation.

Indeed, I could not figure out the role of the two variables “tgv” and “-0.1*c” in the following snippet of code, which I am reporting compared to the theoretical counterpart:

Code line 44 : A = tgv*((lambda + c*(u-pay)) <-0.1*c);

Presentation: Determine A_k :={S : λ_k (S) + c(u_k (S) -φ(S)) < 0}

Moreover, I also failed to understand why, by taking care of the first-order terms using the approximation with the convect operator, we still achieve convergence to the original solution. In the last month, I attempted multiple times to find answers to these questions, but my efforts were in vain. I would be extremely grateful if you could explain the logic behind this, as you are indeed my last hope to comprehend it.

Code:

```
// file BlackScholes2D.edp
int m=40,L=80,LL=80, j=-1, kmax=7;
real T=1, sigmax=0.35, sigmay=0.3, rho=-0.3, r=0.02, K=40, dt=0.01, c=100, tgv=1000;
mesh th=square(m,m,[L*x,LL*y]);
fespace Vh(th,P1);
Vh aux,pay=max(K-max(x,y),0.), u=pay, v,uold,lambda,A=0;
Vh xveloc = -x*r+x*sigmax^2+x*rho*sigmax*sigmay/2,
yveloc = -y*r+y*sigmay^2+y*rho*sigmax*sigmay/2;
varf lambdamat(u,v) = 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) + dx(u)*dy(v))*rho*sigmax*sigmay*x*y/2);
matrix M = lambdamat(Vh,Vh);
varf lambdarhs(unused,v) = int2d(th)( v*convect([xveloc,yveloc],dt,uold)/dt);
problem eq(u,v) = int2d(th)( u*v*(r+1/dt + A)
+ dx(u)*dx(v)*(x*sigmax)^2/2 + dy(u)*dy(v)*(y*sigmay)^2/2
+ (dy(u)*dx(v) + dx(u)*dy(v))*rho*sigmax*sigmay*x*y/2)
- int2d(th)(A*pay*v + v*convect([xveloc,yveloc],dt,uold)/dt)+ on(2,3,u=0);
for (int n=0; n*dt <= T; n++)
{
if(j>T/dt/2) { th = adaptmesh(th,u,verbosity=1,abserror=1,nbjacoby=2,
err=0.0001, 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; v=max(u,pay)-pay; lambda=lambda; pay=pay;
M = lambdamat(Vh,Vh);
plot(v,wait=0,value=1,fill=1);
};
uold=u;
int k=0; aux=0;
while(int2d(th)((u-aux)^2) > 0.8 && k<kmax)
{
aux=u;
eq;
v[] = lambdarhs(0,Vh); lambda[] = M*u[]; lambda[] = v[] -lambda[];
A = tgv*((lambda + c*(u-pay)) <-0.1*c);
k++;
}
cout<<k<<endl; j=j+1;
};
v=max(u,pay)-pay;
plot(v,wait=1,value=1,fill=1);
plot(u,wait=1,value=1,fill=1);
plot(th,wait=1);
```