# Initial condition not used in parablic equation

**URL:** <https://community.freefem.org/t/initial-condition-not-used-in-parablic-equation/2652>\
**Category:** General Discussion\
**Created:** [August 11, 2023, 6:21pm UTC](https://community.freefem.org/t/initial-condition-not-used-in-parablic-equation/2652 "2023-08-11T18:21:18Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Anita](https://avatars.discourse-cdn.com/v4/letter/a/49beb7/32.png) [@Anita](https://community.freefem.org/u/Anita)\
**Post date:** [August 11, 2023, 6:21pm UTC](https://community.freefem.org/t/initial-condition-not-used-in-parablic-equation/2652/1 "2023-08-11T18:21:18Z")

</div>

I’m still struggling with my system of equations but simplified it to 1 parabolic equation to find the issue and feal like I lost some basics.

In below equation value of initial condition doesn’t influence results, looks like only boundary condition changes function values. What’s more function values don’t change much in following time steps.

What I’m missing?

```auto

border C01(t=0.06, 2.08){x=t; y=0.002; }
border C02(t=0.002, 1.002){x=2.08; y=t;}
border C03(t=2.08, 0.06){x=t; y=1.002;}
border C04(t=1.002, 0.002){x=0.06; y=t;}  
int n = 100;
mesh Th = buildmesh(C01(n) + C02(n/2) + C03(n) + C04(n/2));

fespace Vh1(Th,P1);
// Parameters
real Dgamma = 1000;
real beta = 0.85;
real gamma0 = 0.85;
real w1 = 60; //boundary condition

real t=0, T=5, dt=1;
 
Vh1 gamma = gamma0, // initial condition
gammaold,
o2,
uv1 = 1,
uv2 = 1;

macro wektoru(w1,w2) [w1,w2]	//EOM
macro grad(ff) [dx(ff),dy(ff)]	//EOM

problem Probgamma(gamma,o2,solver=LU)=
int2d(Th,qft=qf1pTlump)(gamma/dt*o2)
-int2d(Th,qft=qf1pTlump)(gammaold/dt*o2)
-int2d(Th)(Dgamma*(dx(gamma)*dx(o2) + dy(gamma)*dy(o2))) 
-int2d(Th)(wektoru(uv1, uv2)'*grad(gamma)*o2) 
+int2d(Th)(beta*gamma*o2)

+on(C04, gamma = w1) 
;

for(real t=t;t<=T;t+=dt)
{
cout<<"===================================Solution at time: "<<t<<endl;

gammaold=gamma;
Probgamma; 

plot(gamma,value=1,fill=1,wait=0, cmm="solution gamma at time: "+t);

}

```

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [August 11, 2023, 11:29pm UTC](https://community.freefem.org/t/initial-condition-not-used-in-parablic-equation/2652/2 "2023-08-11T23:29:14Z")

</div>

I didn’t try to write it out but just looking at the problem the laplacian coef is  
large suggesting second der is going to be “small.” It looks like the time derivative goes  
to zero when gamma is just the boundary condition plus an exponential ( almost linear )  
change in the x direction. I played with solver and took out quadrature rule not  
much change although setting “init=1” seems to change things a little… Beta and 1/dt  
seem to be similar magnitude although again you should probably just write it  
out and look at the numbers. I should know the damped thing pretty well but its late here 🙂 There are cases when a condition can be ignored, for example setting  
initial value and initial derivative without a mixed formulation but I’m not  
sure any case like that exists here. What result were your expecting or what equation  
were you trying to implement?

```auto
load "medit"
border C01(t=0.06, 2.08){x=t; y=0.002; }
border C02(t=0.002, 1.002){x=2.08; y=t;}
border C03(t=2.08, 0.06){x=t; y=1.002;}
border C04(t=1.002, 0.002){x=0.06; y=t;}  
int n = 100;
mesh Th = buildmesh(C01(n) + C02(n/2) + C03(n) + C04(n/2));

fespace Vh1(Th,P1);
// Parameters
real Dgamma = 1000;
real beta = 0.85;
real gamma0 = 0.85;
real w1 = 60; //boundary condition

real t=0, T=5*100, dt=1;
 
Vh1 gamma = gamma0, // initial condition
gammaold,
o2,
uv1 = 1,
uv2 = 1;

macro wektoru(w1,w2) [w1,w2]	//EOM
macro grad(ff) [dx(ff),dy(ff)]	//EOM

//problem Probgamma(gamma,o2,solver=LU,init=1)=
problem Probgamma(gamma,o2,solver=sparsesolver,init=0)=
//int2d(Th,qft=qf1pTlump)(gamma/dt*o2)
int2d(Th)(gamma/dt*o2)
//-int2d(Th,qft=qf1pTlump)(gammaold/dt*o2)
-int2d(Th)(gammaold/dt*o2)
-int2d(Th)(Dgamma*(dx(gamma)*dx(o2) + dy(gamma)*dy(o2))) 
-int2d(Th)(wektoru(uv1, uv2)'*grad(gamma)*o2) 
+int2d(Th)(beta*gamma*o2)

+on(C04, gamma = w1) 
;

for(real t=t;t<=T;t+=dt)
{
cout<<"===================================Solution at time: "<<t<<endl;

gammaold=gamma;
Probgamma; 

plot(gamma,value=1,fill=1,wait=0, cmm="solution gamma at time: "+t);
//medit("soln ",Th,gamma);
}

```

---

<div class="post-metadata">

**Author:** ![Anita](https://avatars.discourse-cdn.com/v4/letter/a/49beb7/32.png) [@Anita](https://community.freefem.org/u/Anita)\
**Post date:** [August 12, 2023, 7:08am UTC](https://community.freefem.org/t/initial-condition-not-used-in-parablic-equation/2652/3 "2023-08-12T07:08:06Z")

</div>

Thanks for answering. I’ll look at it deeply again after the weekend.

Checking shortly at your code changes I see that initial condition influences results – good and expected. However, the solution doesn’t change in time like the gamma is constant – so the time derivative is not influencing.

If the initial condition is beta = 0.85 and boundary condition w1 = 60 I would expected values ~between these two in initial step and changing in the time. This equation describes oxygen flow in tissues so expect oxygen from boundary migrates to the whole tissue.

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [August 12, 2023, 10:09am UTC](https://community.freefem.org/t/initial-condition-not-used-in-parablic-equation/2652/4 "2023-08-12T10:09:03Z")

</div>

What eqn are you trying to model?  
AFAICT, as written, for the steady state soln characteristic eqn in x direction roots I got (.5/D)(1 +/- sqrt(1-4BD)).  
Your x values are somewhere around 0-2 or so and D is 1000 so I don’t  
expect the soln to change much but in any case it should not be hard to put in  
analytical soln and compare.

---

<div class="post-metadata">

**Author:** ![Anita](https://avatars.discourse-cdn.com/v4/letter/a/49beb7/32.png) [@Anita](https://community.freefem.org/u/Anita)\
**Post date:** [August 16, 2023, 2:36pm UTC](https://community.freefem.org/t/initial-condition-not-used-in-parablic-equation/2652/5 "2023-08-16T14:36:44Z")

</div>

I’m trying to simulate tumor growth, this is oxygen density equation.  
Probably you gave me important clue which I missed as I’m not so fluent in equations. Original area was rectangle 100-200 not 1-2, I rescaled it due to some other issues I didn’t realize how it influences current version.  
After changing mesh size (to ~100-200) function values start to change however become more unstable in time (growing rapidly or take negative values) – need to verify this deeply.

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [August 16, 2023, 5:23pm UTC](https://community.freefem.org/t/initial-condition-not-used-in-parablic-equation/2652/6 "2023-08-16T17:23:22Z")

</div>

If you are just trying to simulate 1 dimensional diffusion, there are a lot of examples of  
that- heat diffuses and so do charge carriers. Getting a useful result will  
probably require a more complicated model or better boundary conditions.  
Its good to start simple and your current model is so simple there are ways  
to approach it analytically which is always a good sanity check on a computer.  
I’m an engineer so my description of the math may be vulgar but hopefully  
useful 🙂 Probably also you are interested in steady state or quasi static  
situation first where the time derivatives don’t matter so much. But, you  
may want to see analytically what happens to an initial guess of the form  
exp(ikx) or exp(kx) whatever is easier. or for complex k. Then you get  
something like dn/dt=a(k)\*n and you can see if there are conditions where it grows  
indefinitely. If seems to go to zero analytically and the simulation explodes,  
it could be a bug or the time or space scales are too large.