pdepe solver - Time dependent boundary condition
in [Matlab]
|
From: Kamel on 21 Apr 2010 14:25 Hello, I try to impose a boundary condition to resolve a parabolic equation using the pdepe solver. Do you know how to impose a discrete time boundary conition for the right BC like this (series of ramp displacements) : for t<500s, ur = 0.0013*t, 500<t<1000: ur = 0.65, 1000<t<1500: ur = 0.0013*t, 500<t<2000: ur=0.65+0.65. I tried something like: pr = (ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000) + (ur+0.0013*t)*(1000<t<1500)+ (ur+0.65+0.65)*(1500<t<2000); but it works only for the two first blocks (up to t=1000). Please, could you help me to find an easy way to get this BC applied ? Thank you very much, Kamel ps: below is the main program function pdex1 m = 0; x = linspace(0,1.41,50); t = linspace(0,2000,50); sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); % Extract the first solution component as u. u = sol(:,:,1); save -ascii relax_10pts.dat u,x,t % A solution profile can also be illuminating. figure %surf(u,x,t) plot(t,u(:,:)) %plot(0.4*diff(u)) title('Displacement vs time for each z') xlabel('t (s)') ylabel('U (mm)') % -------------------------------------------------------------- function [c,f,s] = pdex1pde(x,t,u,DuDx) c = 1/(0.4*2.7e-3); f = DuDx; s = 0; % -------------------------------------------------------------- function u0 = pdex1ic(x) u0 = 0; % -------------------------------------------------------------- function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)+(ur+0.0013*t)*(1000<t<1500)+(ur+0.65+0.65)*(1500<t<2000); qr = 0;
From: Torsten Hennig on 24 Apr 2010 02:45 > > Thank you Torsten, > > > > dui you know how to use the restart the > calculation > > with the solution obtained so far ? > > > > Thank you very much, > > > > Kamel > > > > An ad-hoc solution: > > Write the solution vector at t=1000 into a > global variable vector and change the value at > the right boundary point appropriately. > Then call pdepe again with different function > handles > for the inital and boundary routines. > In the new routine for the initial values, supply > the > values from the globally accessible solution vector > at t = 1000. > > Best wishes > Torsten. > If this makes technical difficulties, you could modify your original boundary condition for ur to pr = (ur-0.0013*t)*(t<500)+(ur-0.65)*(500<=t<=1000-eps) + (ur-(0.65+0.65/eps*(t-(1000-eps))))* (1000-eps<t<1000)+ + (ur-0.0013*t)*(1000<=t<1500) for a small quantity eps > 0. (Same at t=1500). Best wishes Torsten. > > > > Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> > wrote > > in message > > > <1541508462.1794.1271918912055.JavaMail.root(a)gallium.m > > > athforum.org>... > > > > Hello, > > > > > > > > I try to impose a boundary condition to resolve > a > > > > parabolic equation using the pdepe solver. Do > you > > > > know how to impose a discrete time boundary > > conition > > > > for the right BC like this (series of ramp > > > > displacements) : > > > > > > > > for t<500s, ur = 0.0013*t, 500<t<1000: ur = > 0.65, > > > > 1000<t<1500: ur = 0.0013*t, 500<t<2000: > > > > : ur=0.65+0.65. > > > > > > > > I tried something like: > > > > pr = > (ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000) > > > > > > + (ur+0.0013*t)*(1000<t<1500)+ > > > > (ur+0.65+0.65)*(1500<t<2000); > > > > > > > > > > The boundary condition for u at the right > boundary > > > point does not change continuously at time > t=1000. > > > ur(1000-) = 0.65 ; ur(1000+) = 1.3. > > > Maybe this is a mistake in your setting. > > > If not, you will have to restart pdepe at time > > t=1000 > > > with the solution obtained so far and > artificially > > set > > > u=1.3 in the routine where the initial values > are > > > supplied. > > > > > > Another point is that you have to set > > > pr = (ur-0.0013*t)*(t<500)+(ur-0.65)*(500<t<1000) > > > > + (ur-0.0013*t)*(1000<t<1500)+ > > > (ur-(0.65+0.65))*(1500<t<2000); > > > > > > Best wishes > > > Torsten. > > > > > > > but it works only for the two first blocks (up > to > > > > t=1000). > > > > > > > > Please, could you help me to find an easy way > to > > get > > > > this BC applied ? > > > > > > > > Thank you very much, > > > > > > > > Kamel > > > > ps: below is the main program > > > > > > > > > > > > function pdex1 > > > > > > > > m = 0; > > > > x = linspace(0,1.41,50); > > > > t = linspace(0,2000,50); > > > > > > > > sol = > pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); > > > > % Extract the first solution component as u. > > > > u = sol(:,:,1); > > > > save -ascii relax_10pts.dat u,x,t > > > > > > > > % A solution profile can also be illuminating. > > > > figure > > > > %surf(u,x,t) > > > > plot(t,u(:,:)) > > > > %plot(0.4*diff(u)) > > > > title('Displacement vs time for each z') > > > > xlabel('t (s)') > > > > ylabel('U (mm)') > > > > % > > > > > > > ------------------------------------------------------ > > > > > -------- > > > > function [c,f,s] = pdex1pde(x,t,u,DuDx) > > > > c = 1/(0.4*2.7e-3); > > > > f = DuDx; > > > > s = 0; > > > > % > > > > > > > ------------------------------------------------------ > > > > > -------- > > > > function u0 = pdex1ic(x) > > > > u0 = 0; > > > > % > > > > > > > ------------------------------------------------------ > > > > > -------- > > > > function [pl,ql,pr,qr] = > pdex1bc(xl,ul,xr,ur,t) > > > > pl = ul; > > > > ql = 0; > > > > pr = > > > > > > > ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)+(ur+0.0013 > > > > > *t)*(1000<t<1500)+(ur+0.65+0.65)*(1500<t<2000); > > > > > qr = 0;
From: Kamel on 27 Apr 2010 07:49 Thanks! I will try it. By the way, do you know if I can resolve simultaneously these equations using the pdepe solver: DuDt = k1*k2*DuDx2 DpDx = (1/k1)*DuDt where u, p are dipslacment and pressure fields, k1 and k2 constants. The second equation looks like hyperbolic form and I don't know if we can resolve it with pdepe. Thank you, Kamel Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <989007544.14325.1272105971532.JavaMail.root(a)gallium.mathforum.org>... > > > Thank you Torsten, > > > > > > dui you know how to use the restart the > > calculation > > > with the solution obtained so far ? > > > > > > Thank you very much, > > > > > > Kamel > > > > > > > An ad-hoc solution: > > > > Write the solution vector at t=1000 into a > > global variable vector and change the value at > > the right boundary point appropriately. > > Then call pdepe again with different function > > handles > > for the inital and boundary routines. > > In the new routine for the initial values, supply > > the > > values from the globally accessible solution vector > > at t = 1000. > > > > Best wishes > > Torsten. > > > > If this makes technical difficulties, you could > modify your original boundary condition for ur to > pr = (ur-0.0013*t)*(t<500)+(ur-0.65)*(500<=t<=1000-eps) > + (ur-(0.65+0.65/eps*(t-(1000-eps))))* > (1000-eps<t<1000)+ > + (ur-0.0013*t)*(1000<=t<1500) > for a small quantity eps > 0. > > (Same at t=1500). > > Best wishes > Torsten. > > > > > > > > > > Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> > > wrote > > > in message > > > > > <1541508462.1794.1271918912055.JavaMail.root(a)gallium.m > > > > > athforum.org>... > > > > > Hello, > > > > > > > > > > I try to impose a boundary condition to resolve > > a > > > > > parabolic equation using the pdepe solver. Do > > you > > > > > know how to impose a discrete time boundary > > > conition > > > > > for the right BC like this (series of ramp > > > > > displacements) : > > > > > > > > > > for t<500s, ur = 0.0013*t, 500<t<1000: ur = > > 0.65, > > > > > 1000<t<1500: ur = 0.0013*t, 500<t<2000: > > > > > : ur=0.65+0.65. > > > > > > > > > > I tried something like: > > > > > pr = > > (ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000) > > > > > > > > + (ur+0.0013*t)*(1000<t<1500)+ > > > > > (ur+0.65+0.65)*(1500<t<2000); > > > > > > > > > > > > > The boundary condition for u at the right > > boundary > > > > point does not change continuously at time > > t=1000. > > > > ur(1000-) = 0.65 ; ur(1000+) = 1.3. > > > > Maybe this is a mistake in your setting. > > > > If not, you will have to restart pdepe at time > > > t=1000 > > > > with the solution obtained so far and > > artificially > > > set > > > > u=1.3 in the routine where the initial values > > are > > > > supplied. > > > > > > > > Another point is that you have to set > > > > pr = (ur-0.0013*t)*(t<500)+(ur-0.65)*(500<t<1000) > > > > > > + (ur-0.0013*t)*(1000<t<1500)+ > > > > (ur-(0.65+0.65))*(1500<t<2000); > > > > > > > > Best wishes > > > > Torsten. > > > > > > > > > but it works only for the two first blocks (up > > to > > > > > t=1000). > > > > > > > > > > Please, could you help me to find an easy way > > to > > > get > > > > > this BC applied ? > > > > > > > > > > Thank you very much, > > > > > > > > > > Kamel > > > > > ps: below is the main program > > > > > > > > > > > > > > > function pdex1 > > > > > > > > > > m = 0; > > > > > x = linspace(0,1.41,50); > > > > > t = linspace(0,2000,50); > > > > > > > > > > sol = > > pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); > > > > > % Extract the first solution component as u. > > > > > u = sol(:,:,1); > > > > > save -ascii relax_10pts.dat u,x,t > > > > > > > > > > % A solution profile can also be illuminating. > > > > > figure > > > > > %surf(u,x,t) > > > > > plot(t,u(:,:)) > > > > > %plot(0.4*diff(u)) > > > > > title('Displacement vs time for each z') > > > > > xlabel('t (s)') > > > > > ylabel('U (mm)') > > > > > % > > > > > > > > > > ------------------------------------------------------ > > > > > > > -------- > > > > > function [c,f,s] = pdex1pde(x,t,u,DuDx) > > > > > c = 1/(0.4*2.7e-3); > > > > > f = DuDx; > > > > > s = 0; > > > > > % > > > > > > > > > > ------------------------------------------------------ > > > > > > > -------- > > > > > function u0 = pdex1ic(x) > > > > > u0 = 0; > > > > > % > > > > > > > > > > ------------------------------------------------------ > > > > > > > -------- > > > > > function [pl,ql,pr,qr] = > > pdex1bc(xl,ul,xr,ur,t) > > > > > pl = ul; > > > > > ql = 0; > > > > > pr = > > > > > > > > > > ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)+(ur+0.0013 > > > > > > > *t)*(1000<t<1500)+(ur+0.65+0.65)*(1500<t<2000); > > > > > > > qr = 0;
From: Torsten Hennig on 27 Apr 2010 21:35 > Thanks! I will try it. ... or you can define an event function to restart the simulation at t = 1000. See the pdepe documentation for more details. > By the way, do you know if I can resolve > simultaneously these equations using the pdepe > solver: > DuDt = k1*k2*DuDx2 > DpDx = (1/k1)*DuDt From dp/dx = (du/dt)/k1 = k2*d^2u/dx^2 you get p(x) - p(x_0) = k2 * ((du/dx)(x) - (du/dx)(x_0)) and so p(x) = p(x_0) + k2*((du/dx)(x) - (du/dx)(x_0)) (1) Use pdeval to compute du/dx at the end of the simulation for the output times and 'reconstruct' p by (1). > > where u, p are dipslacment and pressure fields, k1 > and k2 constants. > The second equation looks like hyperbolic form and I > don't know if we can resolve it with pdepe. > > Thank you, > > Kamel Best wishes Torsten.
From: Kamel on 23 May 2010 18:16 Hello Torsten, I would like to prescribe this expression for the right boundary condition (ur vs t) : step_up_times= [100 400 600 800]; n = 0; ur = zeros(1,1007); t = 1:1007; for k = 1:length(ur) if t(k) <= step_up_times(end) && t(k) >= step_up_times(n+1) n = n+1; end ur(k) = .1*n; end This is still a ramp boundary condition. Do you know how to introduce it easily in the function pdex1bc ? function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = "ur vs t" ; qr = 0; Cheers, Kamel Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <947386055.34939.1272432976798.JavaMail.root(a)gallium.mathforum.org>... > > Thanks! I will try it. > > .. or you can define an event function to restart the > simulation at t = 1000. > See the pdepe documentation for more details. > > By the way, do you know if I can resolve > > simultaneously these equations using the pdepe > > solver: > > DuDt = k1*k2*DuDx2 > > DpDx = (1/k1)*DuDt > > From > dp/dx = (du/dt)/k1 = k2*d^2u/dx^2 > you get > p(x) - p(x_0) = k2 * ((du/dx)(x) - (du/dx)(x_0)) > and so > p(x) = p(x_0) + k2*((du/dx)(x) - (du/dx)(x_0)) (1) > > Use pdeval to compute du/dx at the end of the simulation > for the output times and 'reconstruct' p by (1). > > > > where u, p are dipslacment and pressure fields, k1 > > and k2 constants. > > The second equation looks like hyperbolic form and I > > don't know if we can resolve it with pdepe. > > > > Thank you, > > > > Kamel > > Best wishes > Torsten.
|
Pages: 1 Prev: a question about HMM functions in MATLAB Next: Help with Matlab hw! |