| Author | Message | ||
     Geoff Mercer (geoffm) Junior Member Username: geoffm Post Number: 3 Registered: 03-2004 |
Hi, I am not sure if flexpde is able to solve this. I cant seem to get my head around being able to define a localised integral. Here is the problem I have a 1D diffusion system (later to be 2d but lets stick to 1d first) dt(a)=-gamma + D dxx(a) dt(c)=gamma where gamma is a function of a,c and phi. The problem is the phi. phi is defined to be a local integral of c eg phi = integral_(x+delta)^(x-delta) c(z) dz where delta is some given number (not necessarily too small) so how does one specify phi ? Hope you can help. Geoff | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 128 Registered: 06-2003 |
All functions are evaluated within the context of the position where their values is needed - that is, coordinate positions are inserted into the meaning of X, Y, etc according to the evaluation being performed. Definitions with parameters (functions) take an argument list, which supplements current values of coordinates and definitions. So the following should do what you want: Definitions phif(x0) = integral(if abs(x-x0)<delta then c else 0) phi = phif(x) The definition of phi takes the value of the current evaluation coordinate at the point where phi is needed, and passes it to the integrating function phif as x0. The integration sweeps through the mesh cells, passing in X as the current integrand evaluation coordinate, and X0 as the reference position. "c" is evaluated at "x" the current integrand coordinate. This is syntactically correct, at least. I ran a test, but I didn't know what the result was supposed to be, so I'm not certain it did it right. This could get expensive. | ||
| Geoff Mercer (geoffm) Member Username: geoffm Post Number: 4 Registered: 03-2004 |
Robert, that all makes sense but it doesnt seem to work. (see example below). It appears that abs(x-x0) is always = 0 ???? and hence the integral is not localised. In the example below phi should vary with x (since c(x) does) but it is constant across x. When i changed the integral to phif(x0) = bintegral(if abs(x-x0)<delta then x else 0,"the bottom") then it just works out x^2/2]_(0)^(0.5)=1/8 (ie over the entire x domain) instead of the local domain x^2/2]_(x-delta)^(x+delta) = 2*x*delta I hope that you can understand this example. Geoff Mercer {Ants and corpses} title "Ants and Corpse Piles" select smoothinit off errlim 1e-3 cubic { Use Cubic Basis } thermal_colors on variables c a definitions Lx=0.5 Ly=0.05 D=1.3E-3 {Diffusivity} alpha1=31.75 alpha2=1000 alpha3=3.125 alpha4=50 kd=0.75 rho=40 delta=0.01 v=1.6E-2 tend=1 ! phif(x0) = bintegral(if abs(x-x0)<delta then x else 0,"the bottom") phif(x0) = integral(if abs(x-x0)<delta then c else 0) phi = phif(x) gamma=v*( kd*a + alpha1*a*phi/(alpha2+phi) - alpha3*rho*c/(alpha4+phi) ) initial value c=exp(-200*(x-0.25)^2) a=1 equations gamma=dt(c) D*div(grad(a)) -gamma = dt(a) boundaries region 1 START (0,0) line to (Lx,0) line to (Lx,Ly) line to (0,Ly) line to finish start "the bottom" (0,0) line to (Lx,0) finish time 0 to tend by tend/100 monitors for t = 0 by tend/100 to tend elevation(c) from (0,0) to (Lx,0) as "corpes" elevation(a) from (0,0) to (Lx,0) as "ants" elevation(phi) from (0,0) to (Lx,0) as "phi" elevation(gamma) from (0,0) to (Lx,0) as "gamma" histories history(c) at (Lx/2,0) as "c at centre" history(a) at (Lx/2,0) as "a at centre" end | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 130 Registered: 06-2003 |
Unfortunately, it's more complicated than I thought at first. What I said is true in principle, but there is a snag: FlexPDE keeps a scoreboard of valid integrals, to prevent recomputation of expensive things unless their value has changed. This scoreboard is not handled properly in the case where an integral uses the argument of a function call which contains the integral. I dug a little way into it and found that fixing it is really rather messy. So, I'm going to have to defer this effort until a later time, and say that for now, FlexPDE can't do what you want. Sorry. |