Steady states of Nonlinear PDE system

General Discussion
1

Hi,

I have created a model for non-linear PDE model for cancer growth with a time stepping loop. However when setting the initial conditions to the spatially homogenous state, the populations change. This indicates something must be wrong with the code as the populations should remain constant. I believe it may be something to do with my time stepping but I am unsure.

I would appreciate any suggestions on what might be going wrong and how to fix it.

Thanks!

SpaceP1 c, cold, clogvalue, u;
SpaceP1 Ncell, nold, nlogvalue;
SpaceP1 vegfA, aold;
SpaceP1 vegfB, bold;
SpaceP1 X, xold;
SpaceP1 E, eold, Ediv;
SpaceP1 Kc, Kn, Pan, Pac, Pbn, Pbc, f;


varf Cequation(c, v) =
    int2d(Mesh)(c/dt*v)
    + int2d(Mesh)(Dc*(grad(c)' * grad(v)))
    + int2d(Mesh)(cold/dt*v)
    + int2d(Mesh)(-c*log((clogvalue/Kc) + 0.01)*v)
;

varf Nequation(Ncell, vn) =
    int2d(Mesh)(Ncell/dt*vn)
    + int2d(Mesh)(Dn*(grad(Ncell)' * grad(vn)))
    + int2d(Mesh)(nold/dt*vn)
    + int2d(Mesh)(-lambdar*Ncell*log((nlogvalue/Kn) + 0.01)*vn)
;

varf Aequation(vegfA, va) =
    int2d(Mesh)(vegfA/dt*va)
    + int2d(Mesh)(Da*(grad(vegfA)' * grad(va)))
    + int2d(Mesh)(aold/dt*va)
    + int2d(Mesh)(Pan*nold*va)
    + int2d(Mesh)(Pac*cold*va)
    + int2d(Mesh)(-deltaA * vegfA* va)
;

varf Bequation(vegfB, vb) =
    int2d(Mesh)(vegfB/dt*vb)
    + int2d(Mesh)(Db*(grad(vegfB)' * grad(vb)))
    + int2d(Mesh)(bold/dt*vb)
    + int2d(Mesh)(Pbn*nold*vb)
    + int2d(Mesh)(Pbc*cold*vb)
    + int2d(Mesh)(- deltaB * vegfB* vb)
;

varf Xequation(X, vx) =
    int2d(Mesh)(X/dt*vx)
    + int2d(Mesh)(Dx*(grad(X)' * grad(vx)))
    + int2d(Mesh)(xold/dt*vx)
    + int2d(Mesh)(Q*eold*vx)
    + int2d(Mesh)(-Q*X*eold*vx)
    + int2d(Mesh)(-cold*X*vx)
    + int2d(Mesh)(-nold*X*vx)
;

varf Eequation(E, ve) =
    int2d(Mesh)(E/dt*ve)
    + int2d(Mesh)(De*(grad(E)' * grad(ve)))
    + int2d(Mesh)(Cea*(grad(aold)' * grad(ve)))
    + int2d(Mesh)(eold/dt*ve)
    + int2d(Mesh)(f*E*ve)

    + int2d(Mesh)(- deltaE * E* ve)
;

X = 0.58019166;
c = 2.2463892;
Ncell = 0.547827;
vegfA = 10.45252677;
vegfB = 0.70371085;
E = 3.86171686;

Kc = (Kcmin*hc+ Kcmax * X)/(hc + X);
Kn = (Knmin*hn+ Knmax * X)/(hn + X);
Pan = (Panmax * Kan^nan + Panmin* X^nan)/ (X^nan + Kan^nan);
Pac = (Pacmax * Kac^nac + Pacmin* X^nac)/ (X^nac + Kac^nac);
Pbn = (Pbnmin * Kbn^nbn + Pbnmax* X^nbn)/ (X^nbn + Kbn^nbn);
Pbn = (Pbcmin * Kbc^nbc + Pbcmax* X^nbc)/ (X^nbc + Kbc^nbc);
f = (R/(1+exp(vegfB-vegfA+gamma)))*(1-(E/(omega*Ncell+c)));

matrix<real> Ax = Xequation(SpaceP1, SpaceP1, solver=sparsesolver);
real[int] rhsX = Xequation(0, SpaceP1);
matrix<real> A = Cequation(SpaceP1, SpaceP1, solver=sparsesolver);
real[int] rhs = Cequation(0, SpaceP1);
matrix<real> An = Nequation(SpaceP1, SpaceP1, solver=sparsesolver);
real[int] rhsn = Nequation(0, SpaceP1);
matrix<real> Aa = Aequation(SpaceP1, SpaceP1, solver=sparsesolver);
real[int] rhsA = Aequation(0, SpaceP1);
matrix<real> Ab = Bequation(SpaceP1, SpaceP1, solver=sparsesolver);
real[int] rhsB = Bequation(0, SpaceP1);
matrix<real> Ae = Eequation(SpaceP1, SpaceP1, solver=sparsesolver);
real[int] rhsE = Eequation(0, SpaceP1);

for(int i = 0; i < T/dt; i++)
{
  //  cout << "iteration: " << i << endl;

    xold = X;
    cold = c;
    nold = Ncell;
    bold = vegfB;
    aold = vegfA;
    eold = E;

    rhsX = Xequation(0, SpaceP1);
    X[] = Ax^-1 * rhsX;

    Kc = (Kcmin*hc+ Kcmax * xold)/(hc + xold);
	clogvalue = cold;
	rhs  = Cequation(0, SpaceP1);
    c[] = A^-1 * rhs;
    
    Kn = (Knmin*hn+ Knmax * xold)/(hn + xold);
    nlogvalue = nold;
	rhsn  = Nequation(0, SpaceP1);
    Ncell[] = An^-1 * rhsn;

	Pan = (Panmax * Kan^nan + Panmin* xold^nan)/ (xold^nan + Kan^nan);
    Pac = (Pacmax * Kac^nac + Pacmin* xold^nac)/ (xold^nac + Kac^nac);
    rhsA = Aequation(0, SpaceP1);
    vegfA[] = Aa^-1 * rhsA;

    Pbn = (Pbnmin * Kbn^nbn + Pbnmax* xold^nbn)/ (xold^nbn + Kbn^nbn);
    Pbn = (Pbcmin * Kbc^nbc + Pbcmax* xold^nbc)/ (xold^nbc + Kbc^nbc);
    rhsB = Bequation(0, SpaceP1);
    vegfB[] = Ab^-1 * rhsB;
    
    Ediv = eold;
    f = (R/(1+exp(bold-aold+gamma)))*(1-(eold/(omega*nold+cold)));
    rhsE = Eequation(0, SpaceP1);
    E[] = Ae^-1 * rhsE;