My PDE is not obeying conservation principle
in [Matlab]
|
From: Sunipa Som on 1 Jul 2010 04:49 Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <906445592.28720.1277909699309.JavaMail.root(a)gallium.mathforum.org>... > > 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. Hi, r is the radius of a sphere in which particles are moving. Regards, Sunipa som
From: Torsten Hennig on 1 Jul 2010 02:17 > Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote > in message > <906445592.28720.1277909699309.JavaMail.root(a)gallium.m > athforum.org>... > > > 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. > > Hi, > r is the radius of a sphere in which particles are > moving. > > Regards, > Sunipa som Then I think already your principle equation is wrong. Because of the spherical coordinate system, it should read dN/dt + 1/r^2*d/dr(r^2*v_r*N) + F*dN/dp = 0. However: If your equation reads dN/dt + v_x*dN/dx + F/m*dN/dv_x, Wikipedia says that N*dx*dv_x is the number of particles in the volume dx x dv_x. Thus the total number of particles is given by N_sum = int_2d N dx dv_x = int_{x=x_min}^{x=x_max} int_{v_x = v_x_min}^{v_x=v_x_max} N(x,v_x) dx dv_x. In spherical coordinates, it should be given by N_sum = 4*pi*int_{r=r_min}^{r=r_max} int_{v_r=v_r_min}^{v_r=v_r_max} r^2*N(r,v_r) d_v_r d_r. Best wishes Torsten.
From: Torsten Hennig on 1 Jul 2010 02:33 > > Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> > wrote > > in message > > > <906445592.28720.1277909699309.JavaMail.root(a)gallium.m > > > athforum.org>... > > > > 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. > > > > Hi, > > r is the radius of a sphere in which particles are > > moving. > > > > Regards, > > Sunipa som > > Then I think already your principle equation is > wrong. > Because of the spherical coordinate system, > it should read > dN/dt + 1/r^2*d/dr(r^2*v_r*N) + F*dN/dp = 0. > I'm sorry; the equation is indeed dN/dt + v_r*dN/dr + F*dN/dp = 0 as I found in the internet. > However: > If your equation reads > dN/dt + v_x*dN/dx + F/m*dN/dv_x, > Wikipedia says that N*dx*dv_x is the number of > particles in the volume dx x dv_x. > Thus the total number of particles is given by > N_sum = int_2d N dx dv_x = int_{x=x_min}^{x=x_max} > int_{v_x = v_x_min}^{v_x=v_x_max} N(x,v_x) dx dv_x. > > In spherical coordinates, it should be given by > N_sum = 4*pi*int_{r=r_min}^{r=r_max} > int_{v_r=v_r_min}^{v_r=v_r_max} r^2*N(r,v_r) d_v_r > d_r. > > Best wishes > Torsten.
From: Sunipa Som on 1 Jul 2010 07:40 Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <341894495.32310.1277980412575.JavaMail.root(a)gallium.mathforum.org>... > > > Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> > > wrote > > > in message > > > > > <906445592.28720.1277909699309.JavaMail.root(a)gallium.m > > > > > athforum.org>... > > > > > 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. > > > > > > Hi, > > > r is the radius of a sphere in which particles are > > > moving. > > > > > > Regards, > > > Sunipa som > > > > Then I think already your principle equation is > > wrong. > > Because of the spherical coordinate system, > > it should read > > dN/dt + 1/r^2*d/dr(r^2*v_r*N) + F*dN/dp = 0. > > > > I'm sorry; the equation is indeed > dN/dt + v_r*dN/dr + F*dN/dp = 0 > as I found in the internet. > > > > However: > > If your equation reads > > dN/dt + v_x*dN/dx + F/m*dN/dv_x, > > Wikipedia says that N*dx*dv_x is the number of > > particles in the volume dx x dv_x. > > Thus the total number of particles is given by > > N_sum = int_2d N dx dv_x = int_{x=x_min}^{x=x_max} > > int_{v_x = v_x_min}^{v_x=v_x_max} N(x,v_x) dx dv_x. > > > > In spherical coordinates, it should be given by > > N_sum = 4*pi*int_{r=r_min}^{r=r_max} > > int_{v_r=v_r_min}^{v_r=v_r_max} r^2*N(r,v_r) d_v_r > > d_r. > > > > Best wishes > > Torsten. Hi, thank you for your help. Then if I do the integration N_sum=0; for i=1:nr for j=1:nk N_sum=N_sum+(4*pi*(r(i))^2.*N_rad(i,j)); end end Is it right way? With Regards, Sunipa Som
From: Torsten Hennig on 5 Jul 2010 22:37
>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 Start from the equation dN/dt + k*c1*dN/dr-c2*r*dN/dk=0 or dN/dt + div(k*c1*N,-c2*r*N) = 0 to be solved over the rectangle V=[r_min;r_max]x[k_min,k_max]. Integrate over the rectangle V with boundary A to get d/dt int_{V} N dV = int_{V} div(-k*c1*N,c2*r*N) dV. Applying Gauss' integral theorem to the right hand side results in d/dt int_{V} N dV = int_{A} (-k*c1*N,c2*r*N)*n dA where n is the unit normal pointing outwards of the rectangle. So if no particles enter the rectangle over the boundary, the quantity int_{V} N dV = int_{r_min}^{r_max} int_{k_min}^{k_max} N(r,k) dr dk is conserved. This quantity can be approximated from your calculations by sum_{i=1}^{i=N-1} sum_{j=1}^{j=M-1} N(r_i,k_j)* (r_(i+1)-r_(i))*(k_(j+1)-k_(j)). Whether this is the number of molecules in the sphere, you must decide from the physical background of your problem. Best wishes Torsten. |