Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have a PDE function that contains a variable (named ka). I want to plot the result of my PDE in 1 plot using several different ka values. However, I'm not sure about how to write down this equation. So far I can Plot this result but only by using a defined value. Is it possible to change my ka value?

If it is not possible, let say that I dont have a single ka value, how can I insert this value in the middle of the equation?

I put my code below

kw = 10^(-10);(*m/s*)
ka = 10^(-10);(*m/s*)
gw = 10;(*N/m^3*)
m2w = -2*10^(-4);(*coefficient of compressibility of water phase with respect to matric suction*)
m1w = -0.5*10^(-4);(*coefficient of compressibility of water phase with respect to axial stress*)
m2a = 1*10^(-4);(*coefficient of compressibility with air phase with respect to matric suction*)
m1a = -2*10^(-4);(*coefficient of compressibility of air phase with respect to axial stress*)

tmin = 0;(*minimum or time when the pressure is applied*)
tmax = 10^10;(*latest inspected time*)
zmin = 0;(*top elevation*)
zmax = 10;(*bottom elevation*)
uatm = 101;(*kPa-atmospheric air pressure*)
step = 5000;(*The higher the step, the higher the accuracy*)
uwi = 40;(*Initial water pressure*)
uai = 20;(*initial air pressure*)

temp = 293.16;(*temperature in kelvin*)
wa = 28.97*10^(-3);(*kg/mol*)
s = 0.8;
n = 0.5;
r = 8.31432;(*Universal gas constant in j/(mol K*)
g = 10;(*gravity acceleration-m2/s*)
da = ka/g;
cwv = kw/(gw*m2w)
cw = ((1 - m2w/m1w)/(m2w/m1w))
ca = (m2a/m1a)/(1 - (m2a/m1a) - (1 - s)*n/(uatm*m1a))
cav = da/(wa/(r*temp)*m1a*uatm*(1 - m2a/m1a - (1 - s)*n/(uatm*m1a)))
cws = m1w/m2w;
cas = 1/(1 - (m2a/m1a) - n*(1 - s)/(uatm*m1a));
dt = tmax - tmin;
is = 300;(*ImageSize*)


(*Coupled equation*)
eq1 = -Derivative[2, 0][uw][z, t]*cwv - Derivative[0, 1][ua][z, t]*cw + (D[q[t], t]*cws) == Derivative[0, 1][uw][z, t];
eq2 = -Derivative[2, 0][ua][z, t]*cav - Derivative[0, 1][uw][z, t]*ca + (D[q[t], t]*cas) == Derivative[0, 1][ua][z, t];

(*Boundary condition*)
ic1 = {uw[z, tmin] == uwi, Derivative[1, 0][uw][zmin, t] == 0, uw[zmax, t] == 0};
ic2 = {ua[z, tmin] == uai, ua[zmin, t] == 0, Derivative[1, 0][ua][zmax, t] == 0};

sol2 = NDSolve[Flatten[{eq1, eq2, ic1, ic2}], {uw, ua}, {z, zmin, zmax}, {t, tmin,tmax}, MaxSteps -> Infinity, MaxStepSize -> dt/step, InterpolationOrder -> All]

Column[Manipulate[Plot[uw[z, t] /. sol2, {z, zmin, zmax}, PlotRange -> Full, AxesLabel -> {z, uw[z, t]}], {t, tmin, tmax}],Manipulate[LogLinearPlot[(uw[z, t] /. sol2), {t, tmin + 0.1, tmax}, PlotRange -> Full, AxesLabel -> {t, uw[z, t]}], {z, zmin, zmax}]]

Regards Martin Wijaya

share|improve this question

You might try solving with different values of ka and then making an interpolation from that. Alternatively could use ParametricNDSolve but this will be slow due to all the work going on. This can be done as follows. First use your setup but without defining ka. Then one can do:

sol2 = ParametricNDSolve[
  Flatten[{eq1, eq2, ic1, ic2}], {uw, ua}, {z, zmin, zmax}, {t, tmin, 
   tmax}, ka, MaxStepSize -> dt/step, InterpolationOrder -> All]

Now it can be manipulated, sort of. I set relevant parameters so that it will not time out and not try to reevalaute too soon.

 Plot[uw[ka][z, t] /. sol2, {z, zmin, zmax}, PlotRange -> Full, 
  AxesLabel -> {z, uw[z, t]}], {t, tmin, tmax}, {ka, 10^(-11), 
  10^(-9)}, SynchronousInitialization -> False, 
 SynchronousUpdating -> False, ContinuousAction -> False]

Even with these settings it is painfully slow. Possibly there are ways to tune for better performance e.g. by different settings for the NDSolve part.

share|improve this answer
Dear Daniel, How to use ParametricNDSolve?I currently using mathematica 8, and it seems that my mathematica can't recognize this function. – Martin Wijaya Feb 7 '13 at 16:19
Another code that you give somehow also don't give any result (just a blank chart). So, I'm not sure why. Since there is no running sign in the mathematica, I think the process has been completed. – Martin Wijaya Feb 7 '13 at 16:25
ParametricNDSolve is new in version 9. – Daniel Lichtblau Feb 7 '13 at 19:31
You could also plot the sensitivity; that way you get an understanding how a change in ka influences the solution. – user21 Feb 7 '13 at 19:55
How to plot the sensitivity? – Martin Wijaya Feb 8 '13 at 6:12

Your Answer


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

Not the answer you're looking for? Browse other questions tagged or ask your own question.