Prev: replace object properties
Next: approximating the derivative of an array
From: Dave Esc on 14 Jul 2010 07:50
Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <789656615.21906.1279031407787.JavaMail.root(a)gallium.mathforum.org>...
> >
> > Hi again
> >
> > I defined the above terms as follows:
> >
> > global c_inf u k index
> > .
> > .
> > .
> > c_inf = 0;
> > c_inf_array = zeros(100,1);
> > u = 0.0005*ones(3,51,100);
> > index = 0;
> > for k = 1:tstep;
>
> Put index = 0 here (_inside_ instead of outside
> the k-loop). I suspect pdepe already finished
> one time step dt ; so index was 51 and has to be
> reset to 0 for the next time step.
>
> > m = 2;
> > x = 0:.02:1;
> > t = [0 dt/2 dt];
> > sol = pdepe(m,@pde1,@pde1ic,@pde1bc,x,t);
> > u(:,:,k) = sol(:,:,1);
> > c_inf = c_inf - A*dt*(u(3,51,k)-u(3,50,k))/dx;
> > c_inf_array(k) = c_inf;
> > end;
> > %-----------------------------------------------------
> > ---------------
> > function u0 = pde1ic(x)
> > global index c_inf u k
> > index = index + 1;
> > u0 = u(3,index,k);
> > %-----------------------------------------------------
> > ---------------
> >
> > but get error:
> >
> > ??? Attempted to access u(3,52,2); index out of
> > bounds because size(u)=[3,51,100].
> >
> > Error in ==> pde1ic at 4
> > u0 = u(3,index,k);
> > Error in ==> pdepe at 224
> > temp = feval(ic,xmesh(1),varargin{:});
> >
> > Error in ==> pde at 18
> > sol = pdepe(m,@pde1,@pde1ic,@pde1bc,x,t);
> >
> > I dont know how to terminate/set maximum value for
> > _index_ and still keep u0 a scalar!
> > Please help.
> >
> > Many Thanks
> >
>
> Best wishes
> Torsten.

Now works with no errors, however u (and sol) is only showing valid results for the first time step and then repeats the results throughout the loop. I guess u0 still isnt getting updated.

Many thanks
From: Torsten Hennig on 14 Jul 2010 04:24
> Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote
> in message
> <789656615.21906.1279031407787.JavaMail.root(a)gallium.m
> athforum.org>...
> > >
> > > Hi again
> > >
> > > I defined the above terms as follows:
> > >
> > > global c_inf u k index
> > > .
> > > .
> > > .
> > > c_inf = 0;
> > > c_inf_array = zeros(100,1);
> > > u = 0.0005*ones(3,51,100);
> > > index = 0;
> > > for k = 1:tstep;
> >
> > Put index = 0 here (_inside_ instead of outside
> > the k-loop). I suspect pdepe already finished
> > one time step dt ; so index was 51 and has to be
> > reset to 0 for the next time step.
> >
> > > m = 2;
> > > x = 0:.02:1;
> > > t = [0 dt/2 dt];
> > > sol = pdepe(m,@pde1,@pde1ic,@pde1bc,x,t);
> > > u(:,:,k) = sol(:,:,1);
> > > c_inf = c_inf -
> A*dt*(u(3,51,k)-u(3,50,k))/dx;
> > > c_inf_array(k) = c_inf;
> > > end;
> > >
> %-----------------------------------------------------
> > > ---------------
> > > function u0 = pde1ic(x)
> > > global index c_inf u k
> > > index = index + 1;
> > > u0 = u(3,index,k);
> > >
> %-----------------------------------------------------
> > > ---------------
> > >
> > > but get error:
> > >
> > > ??? Attempted to access u(3,52,2); index out of
> > > bounds because size(u)=[3,51,100].
> > >
> > > Error in ==> pde1ic at 4
> > > u0 = u(3,index,k);
> > > Error in ==> pdepe at 224
> > > temp = feval(ic,xmesh(1),varargin{:});
> > >
> > > Error in ==> pde at 18
> > > sol = pdepe(m,@pde1,@pde1ic,@pde1bc,x,t);
> > >
> > > I dont know how to terminate/set maximum value
> for
> > > _index_ and still keep u0 a scalar!
> > > Please help.
> > >
> > > Many Thanks
> > >
> >
> > Best wishes
> > Torsten.
>
> Now works with no errors, however u (and sol) is only
> showing valid results for the first time step and
> then repeats the results throughout the loop. I guess
> u0 still isnt getting updated.
>
> Many thanks

If the initialization of u in the first step
works correctly with the setting u0 = u(3,index,k),
it should also work for the following time steps.
Maybe the diffusion is so slow that you don't see
any changes. Try a bigger value for dt.

Best wishes
Torsten.
From: Torsten Hennig on 14 Jul 2010 22:38
>
> Now works with no errors, however u (and sol) is only
> showing valid results for the first time step and
> then repeats the results throughout the loop. I guess
> u0 still isnt getting updated.
>
> Many thanks

The index k of the time loop was set incorrectly.
With the code below, it should work:

global c_inf u k index
R0 = 0.0015;
S = 2.83E-05;
V = 0.000001;
A = R0*S/V;
dx = 0.02;
tau = 1.7556;
dt = tau/99;
tstep = (tau/dt)+1;
c_inf = 0;
c_inf_array = zeros(101,1);
u = 0.0005*ones(3,51,101);
m = 2;
for k = 1:tstep;
index = 0;
x = 0:.02:1;
t = [0 dt/2 dt];
sol = pdepe(m,@pde1,@pde1ic,@pde1bc,x,t);
u(:,:,k+1) = sol(:,:,1);
c_inf = c_inf - A*dt*(u(3,51,k+1)-u(3,50,k+1))/dx;
c_inf_array(k+1) = c_inf;
end;
%--------------------------------------------------
function u0 = pde1ic(x)
global c_inf u k index
index = index + 1;
u0 = u(3,index,k);
%--------------------------------------------------
function [pl,ql,pr,qr] = pde1bc(xl,ul,xr,ur,t)
global c_inf u k index
betaD = 0.02;
pl = 0;
ql = 1;
pr = betaD*(ur - c_inf);
qr = 1;
%---------------------------------------------------

Best wishes
Torsten.
From: Dave Esc on 22 Jul 2010 13:19
Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <2009041047.32292.1279175918491.JavaMail.root(a)gallium.mathforum.org>...
> >
> > Now works with no errors, however u (and sol) is only
> > showing valid results for the first time step and
> > then repeats the results throughout the loop. I guess
> > u0 still isnt getting updated.
> >
> > Many thanks
>
> The index k of the time loop was set incorrectly.
> With the code below, it should work:
>
> global c_inf u k index
> R0 = 0.0015;
> S = 2.83E-05;
> V = 0.000001;
> A = R0*S/V;
> dx = 0.02;
> tau = 1.7556;
> dt = tau/99;
> tstep = (tau/dt)+1;
> c_inf = 0;
> c_inf_array = zeros(101,1);
> u = 0.0005*ones(3,51,101);
> m = 2;
> for k = 1:tstep;
> index = 0;
> x = 0:.02:1;
> t = [0 dt/2 dt];
> sol = pdepe(m,@pde1,@pde1ic,@pde1bc,x,t);
> u(:,:,k+1) = sol(:,:,1);
> c_inf = c_inf - A*dt*(u(3,51,k+1)-u(3,50,k+1))/dx;
> c_inf_array(k+1) = c_inf;
> end;
> %--------------------------------------------------
> function u0 = pde1ic(x)
> global c_inf u k index
> index = index + 1;
> u0 = u(3,index,k);
> %--------------------------------------------------
> function [pl,ql,pr,qr] = pde1bc(xl,ul,xr,ur,t)
> global c_inf u k index
> betaD = 0.02;
> pl = 0;
> ql = 1;
> pr = betaD*(ur - c_inf);
> qr = 1;
> %---------------------------------------------------
>
> Best wishes
> Torsten.

yes it works now, and matches the analytical solution. The confusion with the boundary condition came about as the analytical solution assumes continuous flux through the sphere into the outside volume which gets rid of the need for a mtc term (betaD); and so the boundary condition used is:
function [pl,ql,pr,qr] = pde2bc(xl,ul,xr,ur,t)
global c_inf
pl = 0;
ql = 1;
pr = ur - c_inf;
qr = 0;
instead. I can finally start altering this to account for other phenomina.
Many Thanks for your help Torsten
First  |  Prev  | 
Pages: 1 2 3
Prev: replace object properties
Next: approximating the derivative of an array