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.

I have a solution that I got from NDSolve which are ua and uw.

eq1=Derivative[2, 0][uw][z, t]*cwv-Derivative[0,1][ua][z,t]*cw == Derivative[0, 1][uw][z, t];
eq2=Derivative[2, 0][ua][z, t]*cav-Derivative[0,1][uw][z,t]*ca == Derivative[0, 1][ua][z, t];

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

result = (ua[z, t] /. sol2) - (atm[z, t]) - (uw[z, t] /. sol2);

I manage to solve this equation and I can plot "result" by using Plot3D.

Plot3D[result, {z, zmin, zmax}, {t, tmin, tmax}, 
       PlotRange -> All, AxesLabel -> {z, t, "result"[z, t]}, 
       BoundaryStyle -> Thick, ImageSize -> 500]

However, I also wants to plot this equation by using manipulate command in Plot but I failed. Here is how I do it

Manipulate[Plot[result, {t, tmin, tmax}, PlotRange -> Full, 
                AxesLabel -> {t, "result"[z, t]}], {z, zmin, zmax}]

or by using this code

Manipulate[result[z,t], {t, tmin, tmax}, PlotRange -> Full, 
           AxesLabel -> {t, "result"[z, t]}], {z, zmin, zmax}]

Can someone please tell me how to plot "result" by using manipulate command?I can do it for each ua[z,t] and uw[z,t]

For the complete set of code

kw = 1*10^(-10);(*m/s*)
ka = 8*10^(-13);(*m/s*)
gw = 10;(*kN/m^3*)
m2w = 2.41*10^(-4);(*coefficient of compressibility of water phase with respect to matric suction*)
m1w = 0.8*10^(-4);(*coefficient of compressibility of water phase with respect to axial stress*)
m2a = 1.11*10^(-4);(*coefficient of compressibility with air phase with respect to matric suction*)
m1a = 0.37*10^(-4);(*coefficient of compressibility of air phase with respect to axial stress*)
p0 = 100;(*applied pressure -kPa*)
tmin = 0;(*minimum or time when the pressure is applied*)
tmax = 5000;(*latest inspected time*)
zmin = 0;(*top elevation*)
zmax = 0.02;(*bottom elevation*)
uatm = 101;(*kPa-atmospheric air pressure*)
step = 1000;(*The higher the step, the higher the accuracy*)
ui = -400;(*Initial matric suction*)

temp = 300;(*temperature in kelvin*)
wa = 28.97;(*kg/kmol*)
s = 0.7887;
n = 1.0696/(1 + 1.0696);
(*CALCULATION*)
r = 8.31432;(*Universal gas constant in j/(mol K*)
g = 9.8;(*gravity acceleration-m2/s*)
da = ka/g;(*function varies with matric suction, havent been defined yet*)
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)*uatm))/(m1a*(1 - m2a/m1a) - (1 - s)*n);
dt = tmax - tmin;

eq1 = Derivative[2, 0][uw][z, t]*cwv - Derivative[0, 1][ua][z, t]*cw == Derivative[0, 1][uw][z, t];
eq2 = Derivative[2, 0][ua][z, t]*cav - Derivative[0, 1][uw][z, t]*ca == Derivative[0, 1][ua][z, t];

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

Column[Plot3D[uw[z, t] /. sol2, {z, zmin, zmax}, {t, tmin, tmax}, PlotRange -> All, AxesLabel -> {z, t, uw[z, t]}, BoundaryStyle -> Thick, ImageSize -> 500], Plot3D[ua[z, t] /. sol2, {z, zmin, zmax}, {t, tmin, tmax}, PlotRange -> All, AxesLabel -> {z, t, ua[z, t]}, BoundaryStyle -> Thick, ImageSize -> 500]]

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

Column[Manipulate[Plot[ua[z, t] /. sol2, {z, zmin, zmax}, PlotRange -> Full, AxesLabel -> {z, ua[z, t]}], {t, tmin, tmax}],Manipulate[Plot[ua[z, t] /. sol2, {t, tmin, tmax}, PlotRange -> Full, AxesLabel -> {t, ua[z, t]}], {z, zmin, zmax}]]

Column[Manipulate[Plot[{uw[z, tmin] /. sol2, uw[z, t1] /. sol2, uw[z, t2] /. sol2, uw[z, t3] /. sol2, uw[z, t4] /. sol2, uw[z, tmax] /. sol2}, {z, zmin, zmax}, AxesLabel -> {z, uw[z, t]}, PlotRange -> Full, AxesLabel -> {z, uw[z, t]}], {t1, tmin + dt/1000, tmin + dt/51}, {t2, tmin + dt/50, tmin + dt/16}, {t3, tmin + dt/15,tmin + dt/3}, {t4, tmin + dt/3, tmin + dt/1}], Manipulate[Plot[{ua[z, tmin] /. sol2, ua[z, t1] /. sol2, ua[z, t2] /. sol2, ua[z, t3] /. sol2, ua[z, t4] /. sol2, ua[z, tmax] /. sol2}, {z, zmin, zmax}, AxesLabel -> {z, ua[z, t]}, PlotRange -> Full], {t1, tmin + dt/500, tmin + dt/51}, {t2, tmin + dt/50, tmin + dt/16}, {t3, tmin + dt/15, tmin + dt/3}, {t4, tmin + dt/3, tmin + dt/1}]]

atm[z, t] = uatm
result == (ua[z, t] /. sol2) - (atm[z, t]) - (uw[z, t] /. sol2);
Plot3D[result, {z, zmin, zmax}, {t, tmin, tmax}, PlotRange -> All, AxesLabel -> {z, t,"result"[z, t]}, BoundaryStyle -> Thick, ImageSize -> 500]
Manipulate[result, {t, tmin, tmax}, PlotRange -> Full, AxesLabel -> {t, "result"[z, t]}], {z, zmin, zmax}]

Regards Martin Wijaya

share|improve this question
add comment

1 Answer

I got your code to work (barring a couple syntax errors) by making result a function of z and t. I copied and pasted everything, but changed the final lines to:

result[z_, t_] = (ua[z, t] /. sol2) - (atm[z, t]) - (uw[z, t] /. sol2);
Plot3D[result[z, t], {z, zmin, zmax}, {t, tmin, tmax}, PlotRange -> All, AxesLabel -> {z, t, "result"[z, t]}, BoundaryStyle -> Thick, ImageSize -> 500]
Manipulate[Plot[result[z, t], {t, tmin, tmax}, PlotRange -> Full, AxesLabel -> {t, "result"[z, t]}], {z, zmin, zmax}]
share|improve this answer
    
Wow...that is really nice...thanks –  Martin Wijaya Feb 6 '13 at 3:18
    
@MartinWijaya The reason why this works and your original did not is that Manipulate localizes variables by default, which means the local z in the Manipulate is a different variable than the global z in result. This solution gets around that by passing the local value in the function call. –  Michael E2 Feb 6 '13 at 15:59
add comment

Your Answer

 
discard

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.