Insight on FreeFem++ Script for American Basket

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);

Hi, not wanting to start the new year with doubt, I have rewritten along the lines of the VI-adap.edp example the variational inequality of the basket price of American options. Given my inexperience, would anyone be able to tell me if this is a good idea?
The results obtained are similar to the code found on the internet which had a slightly different implementation from the theory. Could you please check to see if I have written any corruptions?
Thank you and happy new year.

//--- algo parameters
int m=40;
int L=80;
int LL=80;
int j=-1;
int kmax=1;
int k=0;

//--- financial parameters
real T=1;
real sigmax=0.35;
real sigmay=0.3;
real rho=-0.3;
real r=0.02;
real K=40;
real dt=0.01;
//--- 

bool debug = false;  
//mesh Th=square(20,20);
 mesh Th=square(m,m,[L*x,LL*y]);
 real eps=0.35;
 fespace Vh(Th,P1);         // P1 FE space
 int n = Vh.ndof;           // number of Degree of freedom
 Vh uh,uhp,vh;              // solution and previous one
 Vh Ik;                     // to def the set where the containt is reached. 
 real[int] rhs(n);          // to store the right and side of the equation 
 real c=1000;               // the parameter of the algoritme
 func f=1;                  // right hand side function 
 func fd=0;                 // Dirichlet   boundary condition function
 Vh g=1;
// array to store   
real[int] Aii(n),Aiin(n);    // store the diagonal of the matrix

real tol=0.001,tolmin=0.0001;
real tgv = 1e30;              // a huge value of exact penalisation of boundary condition
real res=0.;
// American option early exercise and other 
func e = max(K-max(x,y),0.);
Vh pay=e;
uh=pay;
Vh xveloc = -x*r+x*sigmax^2+x*rho*sigmax*sigmay/2;
Vh yveloc = -y*r+y*sigmay^2+y*rho*sigmax*sigmay/2;

//  the variational form of the problem:
varf a(uh,vh) = int2d(Th)(  uh*vh*(r+1/dt)
    + dx(uh)*dx(vh)*(x*sigmax)^2/2 + dy(uh)*dy(vh)*(y*sigmay)^2/2
    + (dy(uh)*dx(vh) + dx(uh)*dy(vh))*rho*sigmax*sigmay*x*y/2)
    - int2d(Th)(vh*convect([xveloc,yveloc],dt,uhp)/dt)
    + on(2,3,uh=0);


// two version of the problem  
matrix A=a(Vh,Vh,tgv=tgv,solver=CG);
matrix AA=a(Vh,Vh);

//  the mass Matrix construction: 
varf vM(uh,vh) = int2d(Th)(uh*vh);
matrix M=vM(Vh,Vh); // to do a fast computing of $L^2$ norm : sqrt( u'*(w=M*u)) 
Aii=A.diag; // get the diagonal of the matrix 

rhs = a(0,Vh,tgv=tgv);
Ik =0;
uhp=0;
Vh lambda=0;
int kadapt=0,kkadapt=0, iter=0;
real[int] b(n) ;
real[int] Ak(n);            //  the complementary of Ik ( !Ik = (Ik-1))
j=0;

// iterate for all time T
// Start counting time
real time=clock();
while (iter*dt <= T)
{ 
  k=0;
  while(k<kmax)
  {
    rhs = a(0,Vh,tgv=tgv);
    b=rhs;     //  get a copy of the Right hand side 
    // The operator Ik- 1. is not implement so we do:
    Ak= 1.; Ak  -= Ik[];        // build Ak  = ! Ik 
   
    //  adding new locking  condition on b and on the diagonal if (Ik ==1 )
    b = Ik[] .*pay[];
    b *= tgv;
    b  -=  Ak .* rhs;
    Aiin = Ik[] *  tgv;
    Aiin  +=  Ak  .* Aii;       //set  Aii= tgv  $ i \in Ik $
    A.diag = Aiin;              //  set the matrix diagonal 
    set(A,solver=CG);           // important to change precondiconning  for solving
    
    //SOLVE THE PROBLEM
    uh[] = A^-1* b;             //  solve the problem with more locking condition
  
    //plot(uh, cmm = "uh after solver", wait=debug);
    
    lambda[] = AA * uh[];       //  compute the residual ( fast with matrix)
    lambda[] = -rhs-lambda[];   // remark rhs -

    Ik = (lambda + c*(uh-pay)) < 0;  //  set the new value 

    //plot(Ik, wait=debug,cmm=" lock set ",value=1 );
    //plot(uh,wait=debug,cmm="uh");
    
    // trick to compute  $L^2$ norm of the variation
    real[int] diff(n),Mdiff(n);  
    diff= uh[]-uhp[];    
    Mdiff = M*diff; 
    real err = sqrt(Mdiff'*diff);
    cout << " ----- err norm L2 " << err << " ----- kkadapt = " << kkadapt <<endl;
    if(err<eps && kkadapt ) 
    {
        break;              //stop test
       cout <<" STOP test at k = "<<k<<endl;
    }
    bool adapt = err<eps || (j>T/dt/4); //sure for each ...
    if(adapt && iter>1)
    { 
      kadapt++;
      j=0;
      Th=adaptmesh(Th,uh,err=tol);
      kkadapt = tol == tolmin; // we reacht  the bound       
      tol = max(tol/2,tolmin);
      cout << " ++ tol = " << tol << " kadapt = " << kadapt << " reachted bound? " << kkadapt <<endl;
      //plot(Th,wait=debug);
      //plot(uh, wait=debug);
      Vh xveloc = -x*r+x*sigmax^2+x*rho*sigmax*sigmay/2;
      Vh yveloc = -y*r+y*sigmay^2+y*rho*sigmax*sigmay/2;
      vh=max(uh,pay)-pay;
      //plot(vh, wait=debug, cmm= "vh");
      uhp=uhp;
      n=uhp.n;	
      uh=uh;
      Ik=Ik;
      lambda=lambda;
      pay=e;
      A=a(Vh,Vh,tgv=tgv,solver=CG);
      AA=a(Vh,Vh);
      M=vM(Vh,Vh);
      Aii.resize(n);	
      Aiin.resize(n);
      Ak.resize(n);
      b.resize(n);
      rhs.resize(n); 
      Aii=A.diag; // get the diagonal of the matrix 
    }
    k++;
    cout<<"k = "<<k<<" iter = "<<iter<<endl;
  }
  j++;
  iter++; 
  cout << "iter = "<< iter <<endl; 
  uhp[]=uh[] ; // set the previous solution 
} 
cout << endl;
cout << "===============================================" << endl;
cout << "====            CPU time                  =====" << endl;
cout << "===============================================" << endl;
cout << " ALL solving steps :::: "  << clock()-time << endl;  
cout<<endl;
cout << "====          uh(0,0)"<<uh(0,0)<<"           =====" << endl;
cout<<endl;
cout << "====          uh(35,41)"<<uh(35,41)<<"       =====" << endl;
vh=max(uh,pay)-pay;
plot(vh,wait=1,value=1,fill=1,cmm="vh");
plot(vh,bb=[[20,20],[50,50]],value=1, fill=0,cmm="vh-zoomed",wait=1,nbiso=50);
plot(uh,wait=1,value=1,fill=1,cmm="uh");
plot(uh,wait=1,value=1,fill=0,cmm="uh", nbiso=50);
plot(Th,cmm="Th",wait=1);
plot(Th,bb=[[20,20],[50,50]],cmm="Th-zoomed",wait=1);
savemesh(Th,"mm",[x,y,uh*10]);
//  compute a cut 
int i;
real[int] xx(L),yy(LL);
for (i=0;i<L;i++)
 {
   x=i/L; y=i/LL;
   xx[i]=i;
   yy[i]=uh(i,i); // value of uh at point (i/L , i/LL) 
 }
 plot([xx,yy],bb=[[0,0],[80,80]], cmm="u like 1d", value = true, ps="VIA-american.eps",wait=true, aspectratio=true); //  like gnuplot plot a cut of u
{ // File for gnuplot
    ofstream gnu("plotVIA-american.gp");
    for (int i = 0; i < L; i++)
        gnu << xx[i] << " " << yy[i] << endl;
}

if you search the documentation for “tgv” you can find the details there are a lot
but not too many instances. This is something like ( pardon my french lol )
tres grande value ( lol ) or the value or penalty used to enforce value boundary
conditions. There are special values to zero out the equation or IIRC do more
complete exact handling.

Thank you very much Mike,
the problem lies in the fact that the tgv is set to 1000 and not as usual to 1e30. Also in the script the theoretical plant is not respected, but the final result seems to be correct.
For days I tried to find an explanation and in the end fearing that I was facing a very clever trick I rewrote the variational inequality of the American option price with the help of the VI-adap.edp example. My fear, being a novice with freefem is that I have written mathematical nonsense. For example, I have already noticed a small error in the code I attached.

Andrea.

I guess I would start with small dimensions and dump or display
the intermediate results that you should be able to check bv hand.
There are cases where ff uses some simple things like “x” as a global
variable and it can create some confusion among other issues.

It did look like you had unreachable code after a break in one of the
posts.

Not sure why you would use a tgv around 1000 but again you can dump the matrix
and rhs and see if it makes sense.

Concerning the value of tgv, there is usually no straightforward answer which value of tgv is best. If you want to look for some literature on this, check out e.g. EUDML  |  Penalties, Lagrange multipliers and Nitsche mortaring where s = 1 / tgv, i.e. look for Nitsche penalization, which includes constraints in the variational formation by adding a term to penalize the violation of the constraint. From a mathematical point of view, you should increase the the penalization parameter and let it go to infinity to be able to ensure that asymptotically the constraint is enforced. From a computational point of view, using a penalty parameter (tgv) which is too large could also lead to ill-posed systems which are harder to solve and in the worst case you could have a solution that only satisfies the constraint but not the original problem itself. (You can find a short discussion of tgv = \gamma on page 1372 in https://link.springer.com/content/pdf/10.1007/s10543-021-00854-3.pdf) Therefore, in my opinion you could play around with different values for tgv, e.g. 0, 10, 100, 1000, 10000, to get a feeling what the solution looks like and how much the constraint is satisfied/violated and then choose a value of tgv that is large enough. I think that for many applications a tgv between 10 and 1000 should be sufficient.