Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

This is potentially a daft question, but I thought I'd ask it; I have some material free to diffuse in a boundary between rn and ro; I've been able to get it working nicely for neumann type boundary conditions, and the relevant section of code is here; q[r] is a defined funtion, and all other constants specified and it all works rather nicely;

eqn = D[Ef[r, t], t] - Def5*( D[Ef[r, t], r, r] + (2/(r )) D[Ef[r, t], r])
 +  q[r ]Ef[r, t];


ic = {Ef[r, 0] == 100, Derivative[1, 0][Ef][rn, t] == 0, 
Derivative[1, 0][Ef][ro, t] == 0};
s = NDSolve[{eqn == 0}~Join~ic, Ef, {r, rn, ro}, {t, 0, 60}];

But here's my question; what if the material was free to diffuse out of the ro boundary into the infinite medium around it? The confinement at rn still holds of course; so if this were the case, Ef[Infinity,t] ==0. Mathematica doesn't like this, so I'm using a real valued high value to simulate such behaviour. My initial conditions now look like

ic = {Ef[r, 0] == 100, Derivative[1, 0][Ef][rn, t] == 0,  Ef[ro + 100000, t] == 0};

This produces a warning that boundary and initial conditions are inconsistent, which of course they are; Ef[r,0] == 100 only holds between the boundaries of ro and rn, and is initially zero beyond ro. Is there any way to put this condition into NDSolve?

Essentially, I want to tell Mathematica that the initial conditions ONLY applies between rn and ro, and that initially, Ef[r > ro,t] == 0 ; I should add I'll still only be plotting the solutions over rn-ro, but I imagine this should fall over time as material diffuses beyond ro and I'd like to see that behaviour!

Thanks!

Edit--- The initial conditions and functions are

Def5 = 5.5*10^-11; Do2 = 2*10^-9;  f = 2.6835*10^-7;
po = 100; ro = 500*10^-6;
a = 5.9336*10^-6; rl = Sqrt[(6*f*Do2*po)/a];
rn  = ro*(0.5 - Cos[(ArcCos[1 - (2*(rl^2))/(ro^2)] - 2*Pi)/3]);
const1 = a/(6*f*Do2); const2 = 2*(rn^3);
p[r_] := po + const1 (r^2 - ro^2 + const2 (1/r - 1/ro));

kd = 0.049; kmn = 1.9 ; kme = 3.45; qm = 9.34*10^-4;

q[r_] := qm (((p[r])/(p[r] + kmn)) ((kme)/(kme + p[r])) + kd);
share|improve this question
    
Can you give examples of numeric values for Def5, rn, ro and q[t] so that we can try ? –  andre Apr 25 '13 at 12:29
    
Sure; all the functions are; Def5 = 5.5*10^-11; Do2 = 2*10^-9; f = 2.6835*10^-7; po = 100; ro = 500*10^-6; a = 5.9336*10^-6; rl = Sqrt[(6*f*Do2*po)/a]; rn = ro*(0.5 - Cos[(ArcCos[1 - (2*(rl^2))/(ro^2)] - 2*Pi)/3]); const1 = a/(6*f*Do2); const2 = 2*(rn^3); p[r_] := po + const1 (r^2 - ro^2 + const2 (1/r - 1/ro)); kd = 0.049; kmn = 1.9 ; kme = 3.45; qm = 9.34*10^-4; q[r_] := qm (((p[r])/(p[r] + kmn)) ((kme)/(kme + p[r])) + kd); –  DRG Apr 25 '13 at 12:37
    
^The above wasn't clear so added it to original post! –  DRG Apr 25 '13 at 12:39
    
I have trouble to make your example work. Ef stays in the interval 100 +/- 10E-13. –  andre Apr 25 '13 at 13:03
    
Is it a problem if we do the approximation that q[r]==0 for r>ro ? –  andre Apr 25 '13 at 13:12
show 1 more comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.