My PDE is not obeying conservation principle
in [Matlab]
|
From: Torsten Hennig on 29 Jun 2010 02:40 > Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote > in message > <1225209501.16645.1277733046680.JavaMail.root(a)gallium. > mathforum.org>... > > > > Hi, > > > > My PDE is dN/dt=-c1*dN/dr+c2*dN/dk > > > > from that equation I am getting solution N(r,k) > for > > > a > > > > particular t. Then I tried to calculate the > total > > > > number of particles for different t by > > > > N_sum=dbl integration over > > > 4*pi*r^2*4*pi*k^2*N(r,k) > > > > dr dk > > > > But it is not coming constant. From time steps > 1 to > > > 5 > > > > it is decreasing and then it is constant. > > > > Can you tell me the procedure of calculating > total > > > > number of particles is correct? or any other > idea? > > > > Thank you > > > > S Som > > > > > > Did you take into account the flux of particles > > > over the boundary of your domain of integration > > > [r;R] x [k,K] ? > > > N will only be constant if the total flux sums to > > > zero. > > > > > > > Of course if meant here: > > N_sum will only be constant if the total flux over > the > > boundaries sums to zero. > > > > > Best wishes > > > Torsten. > > Hi, > Thank you for your help. Now how do I calculate total > flux for my case? > And another thing to calculate > N_sum=dbl integration over 4*pi*r^2*4*pi*k^2*N(r,k) > dr dk > I am doing in this way below > > N_sum=0; > for i=1:nr > for j=1:nk > > > N_sum=N_sum+((r(i))^2.*(k(j))^2.*N_rad(i,j)); > end > end > Is it right? > > With Regards, > Sunipa Som Your PDE reads dN/dt = -c1*dN/dr + c2*dN/dk. Assuming c1 and c2 are constant and do not depend on N, r or k, integrating both sides over the domain V = [r_min,r_max] x [k_min,k_max] gives d/dt int_{V} N dV = int_{V} (-c1*dN/dr + c2*dN/dk) dV = int_{A} (-c1*N,c2*N)^t * n dA where ^t denotes "transposed", n is the normal pointing outwards your domain of integration and A is the boundary of V. The last equality uses Gauss' integral theorem. int_{A} (-c1*N,c2*N)^t * n dA is the flux of N over the boundary of V. If this surface integral is not zero, you cannot expect that the quantity N is conserved. Best wishes Torsten.
From: Torsten Hennig on 29 Jun 2010 03:05 > > Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> > wrote > > in message > > > <1225209501.16645.1277733046680.JavaMail.root(a)gallium. > > > mathforum.org>... > > > > > Hi, > > > > > My PDE is dN/dt=-c1*dN/dr+c2*dN/dk > > > > > from that equation I am getting solution > N(r,k) > > for > > > > a > > > > > particular t. Then I tried to calculate the > > total > > > > > number of particles for different t by > > > > > N_sum=dbl integration over > > > > 4*pi*r^2*4*pi*k^2*N(r,k) > > > > > dr dk > > > > > But it is not coming constant. From time > steps > > 1 to > > > > 5 > > > > > it is decreasing and then it is constant. > > > > > Can you tell me the procedure of calculating > > total > > > > > number of particles is correct? or any other > > idea? > > > > > Thank you > > > > > S Som > > > > > > > > Did you take into account the flux of particles > > > > > over the boundary of your domain of > integration > > > > [r;R] x [k,K] ? > > > > N will only be constant if the total flux sums > to > > > > zero. > > > > > > > > > > Of course if meant here: > > > N_sum will only be constant if the total flux > over > > the > > > boundaries sums to zero. > > > > > > > Best wishes > > > > Torsten. > > > > Hi, > > Thank you for your help. Now how do I calculate > total > > flux for my case? > > And another thing to calculate > > N_sum=dbl integration over > 4*pi*r^2*4*pi*k^2*N(r,k) > > dr dk > > I am doing in this way below > > > > N_sum=0; > > for i=1:nr > > for j=1:nk > > > > > > > > N_sum=N_sum+((r(i))^2.*(k(j))^2.*N_rad(i,j)); > > end > > end > > Is it right? > > > > With Regards, > > Sunipa Som > > Your PDE reads > dN/dt = -c1*dN/dr + c2*dN/dk. > Assuming c1 and c2 are constant and do not depend > on N, r or k, integrating both sides over the > domain V = [r_min,r_max] x [k_min,k_max] gives > d/dt int_{V} N dV = > int_{V} (-c1*dN/dr + c2*dN/dk) dV = > int_{A} (-c1*N,c2*N)^t * n dA > where ^t denotes "transposed", n is the normal > pointing outwards your domain of integration and A is > the > boundary of V. > The last equality uses Gauss' integral theorem. > > int_{A} (-c1*N,c2*N)^t * n dA > is the flux of N over the boundary of V. > > If this surface integral is not zero, you cannot > expect that the quantity N is conserved. > Better: If this surface integral is not zero, you can not expect int_{V} N dV to remain constant over time. > Best wishes > Torsten.
From: Torsten Hennig on 29 Jun 2010 22:32 > Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote > in message > <1225209501.16645.1277733046680.JavaMail.root(a)gallium. > mathforum.org>... > > > > Hi, > > > > My PDE is dN/dt=-c1*dN/dr+c2*dN/dk > > > > from that equation I am getting solution N(r,k) > for > > > a > > > > particular t. Then I tried to calculate the > total > > > > number of particles for different t by > > > > N_sum=dbl integration over > > > 4*pi*r^2*4*pi*k^2*N(r,k) > > > > dr dk > > > > But it is not coming constant. From time steps > 1 to > > > 5 > > > > it is decreasing and then it is constant. > > > > Can you tell me the procedure of calculating > total > > > > number of particles is correct? or any other > idea? > > > > Thank you > > > > S Som > > > > > > Did you take into account the flux of particles > > > over the boundary of your domain of integration > > > [r;R] x [k,K] ? > > > N will only be constant if the total flux sums to > > > zero. > > > > > > > Of course if meant here: > > N_sum will only be constant if the total flux over > the > > boundaries sums to zero. > > > > > Best wishes > > > Torsten. > > Hi, > Thank you for your help. Now how do I calculate total > flux for my case? > And another thing to calculate > N_sum=dbl integration over 4*pi*r^2*4*pi*k^2*N(r,k) > dr dk One question: You seem to work in a double-spherical coordinate system, but your PDE is not written in spherical coordinates. Why is this so ? Usually one would expect your PDE to be dN/dt + 1/r^2 * d/dr(r^2*u1*N) + 1/k^2 * d/dk(k^2*u2*N)=0 where u1, u2 are the velocities in r and k-direction. > I am doing in this way below > > N_sum=0; > for i=1:nr > for j=1:nk > > > N_sum=N_sum+((r(i))^2.*(k(j))^2.*N_rad(i,j)); > end > end > Is it right? > > With Regards, > Sunipa Som Best wishes Torsten.
From: Sunipa Som on 30 Jun 2010 04:15 Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <1178231572.26629.1277879565094.JavaMail.root(a)gallium.mathforum.org>... > > Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote > > in message > > <1225209501.16645.1277733046680.JavaMail.root(a)gallium. > > mathforum.org>... > > > > > Hi, > > > > > My PDE is dN/dt=-c1*dN/dr+c2*dN/dk > > > > > from that equation I am getting solution N(r,k) > > for > > > > a > > > > > particular t. Then I tried to calculate the > > total > > > > > number of particles for different t by > > > > > N_sum=dbl integration over > > > > 4*pi*r^2*4*pi*k^2*N(r,k) > > > > > dr dk > > > > > But it is not coming constant. From time steps > > 1 to > > > > 5 > > > > > it is decreasing and then it is constant. > > > > > Can you tell me the procedure of calculating > > total > > > > > number of particles is correct? or any other > > idea? > > > > > Thank you > > > > > S Som > > > > > > > > Did you take into account the flux of particles > > > > over the boundary of your domain of integration > > > > [r;R] x [k,K] ? > > > > N will only be constant if the total flux sums to > > > > zero. > > > > > > > > > > Of course if meant here: > > > N_sum will only be constant if the total flux over > > the > > > boundaries sums to zero. > > > > > > > Best wishes > > > > Torsten. > > > > Hi, > > Thank you for your help. Now how do I calculate total > > flux for my case? > > And another thing to calculate > > N_sum=dbl integration over 4*pi*r^2*4*pi*k^2*N(r,k) > > dr dk > > One question: You seem to work in a double-spherical > coordinate system, but your PDE is not > written in spherical coordinates. > Why is this so ? > Usually one would expect your PDE to be > dN/dt + 1/r^2 * d/dr(r^2*u1*N) + 1/k^2 * d/dk(k^2*u2*N)=0 > where u1, u2 are the velocities in r and k-direction. > > > I am doing in this way below > > > > N_sum=0; > > for i=1:nr > > for j=1:nk > > > > > > N_sum=N_sum+((r(i))^2.*(k(j))^2.*N_rad(i,j)); > > end > > end > > Is it right? > > > > With Regards, > > Sunipa Som > > Best wishes > Torsten. Hi, My equation is Boltzman equation. Without collision term it is dN/dt+v*dN/dr+F*dN/dp=0 where v is velocity, F is force and p is momentum by rearranging few terms we are getting this dN/dt + k*c1*dN/dr-c2*r*dN/dk=0 here unit of k is 1/meter and r is also length. so, then my procedure of calculating total number of particles is right or I have to do in other way? With Regards, Sunipa Som
From: Torsten Hennig on 30 Jun 2010 06:54
> Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote > in message > <1178231572.26629.1277879565094.JavaMail.root(a)gallium. > mathforum.org>... > > > Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> > wrote > > > in message > > > > <1225209501.16645.1277733046680.JavaMail.root(a)gallium. > > > mathforum.org>... > > > > > > Hi, > > > > > > My PDE is dN/dt=-c1*dN/dr+c2*dN/dk > > > > > > from that equation I am getting solution > N(r,k) > > > for > > > > > a > > > > > > particular t. Then I tried to calculate the > > > total > > > > > > number of particles for different t by > > > > > > N_sum=dbl integration over > > > > > 4*pi*r^2*4*pi*k^2*N(r,k) > > > > > > dr dk > > > > > > But it is not coming constant. From time > steps > > > 1 to > > > > > 5 > > > > > > it is decreasing and then it is constant. > > > > > > Can you tell me the procedure of > calculating > > > total > > > > > > number of particles is correct? or any > other > > > idea? > > > > > > Thank you > > > > > > S Som > > > > > > > > > > Did you take into account the flux of > particles > > > > > over the boundary of your domain of > integration > > > > > [r;R] x [k,K] ? > > > > > N will only be constant if the total flux > sums to > > > > > zero. > > > > > > > > > > > > > Of course if meant here: > > > > N_sum will only be constant if the total flux > over > > > the > > > > boundaries sums to zero. > > > > > > > > > Best wishes > > > > > Torsten. > > > > > > Hi, > > > Thank you for your help. Now how do I calculate > total > > > flux for my case? > > > And another thing to calculate > > > N_sum=dbl integration over > 4*pi*r^2*4*pi*k^2*N(r,k) > > > dr dk > > > > One question: You seem to work in a > double-spherical > > coordinate system, but your PDE is not > > written in spherical coordinates. > > Why is this so ? > > Usually one would expect your PDE to be > > dN/dt + 1/r^2 * d/dr(r^2*u1*N) + 1/k^2 * > d/dk(k^2*u2*N)=0 > > where u1, u2 are the velocities in r and > k-direction. > > > > > I am doing in this way below > > > > > > N_sum=0; > > > for i=1:nr > > > for j=1:nk > > > > > > > > > > N_sum=N_sum+((r(i))^2.*(k(j))^2.*N_rad(i,j)); > > > end > > > end > > > Is it right? > > > > > > With Regards, > > > Sunipa Som > > > > Best wishes > > Torsten. > > Hi, > My equation is Boltzman equation. Without collision > term it is > dN/dt+v*dN/dr+F*dN/dp=0 > where v is velocity, F is force and p is momentum > by rearranging few terms we are getting this > dN/dt + k*c1*dN/dr-c2*r*dN/dk=0 > here unit of k is 1/meter and r is also length. > so, then my procedure of calculating total number of > particles is right or I have to do in other way? > > With Regards, > Sunipa Som Is r just a length coordinate ? Or should it be the radius of a sphere in which particles are moving ? Best wishes Torsten. |