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 2 variables as a results from the NDSolve[] which are ua and uw. This is how I do it, and I got the result for ua and uw just fine and I can plot them.

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];`

However, I want to subtract the two variables such that: result = (ua-101)-uw and I failed. I use this as a code to subtract

result = ua[z, t] /. sol2 - 101 - uw[z, t] /. sol2;

Can someone please tell me the reason why it is failed and what should I do to fix this?

The complete code is given as follows:

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 = 200;(*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 = -300;(*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;

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;

(*Assuming that dua/dt is insignificant, Consolidation equation for unsaturated soil is duw/dt=cwv(d^2uw/dy^2)*)

eq = Derivative[2, 0][uw][z, t]*cwv == Derivative[0, 1][uw][z, t];

ic = {uw[zmax, t] == ui + p0, uw[z, tmin] == ui};

sol = NDSolve[Flatten[{eq, ic}], uw, {z, zmin, zmax}, {t, tmin, tmax},MaxSteps -> Infinity, MaxStepSize -> dt/step,InterpolationOrder -> All]

is = 300;(*ImageSize*)

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

Column[Manipulate[Plot[uw[z, t] /. sol, {z, zmin, zmax}, PlotRange -> All, AxesLabel -> {z, uw[z, t]}, ImageSize -> is], {t, tmin, tmax}], Manipulate[Plot[uw[z, t] /. sol, {t, tmin, tmax}, PlotRange -> All, AxesLabel -> {t, uw[z, t]}, ImageSize -> is], {z, zmin, zmax}],  Manipulate[Plot[{uw[z, tmin] /. sol, uw[z, t1] /. sol, uw[z, t2] /. sol, uw[z, t3] /. sol, uw[z, t4] /. sol, uw[z, tmax] /. sol}, {z, zmin,zmax}, AxesLabel -> {z, uw[z, t]}, ImageSize -> is], {t1, tmin + dt/300, tmin + dt/51}, {t2, tmin + dt/50, tmin + dt/16}, {t3, tmin + dt/15, tmin + dt/3}, {t4, tmin + dt/3, tmin + dt/1}]]


(*Assuming that dua/dt is significant*)
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];


(*Coupled equation version*)
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}]]

(*Comparison*)

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

(*Matric Suction*)

result = ua[z, t] /. sol2 - 101 - uw[z, t] /. sol2;
Plot3D[result, {z, zmin, zmax}, {t, tmin, tmax}, PlotRange -> All, AxesLabel -> {z, t, \[Psi][z, t]}, BoundaryStyle -> Thick, ImageSize -> 500]    `

Regards,

share|improve this question
    
This question might be useful for future visitors if you strip it down to essentials (in fact the bit @AlbertRetey points out for you). Please consider editing. –  Yves Klett Feb 8 '13 at 7:50
add comment

1 Answer

I think it's just a matter of operator precedence. It should work when you add parenthesis at the right places:

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

as ReplaceAll (== /.) will apply the replacement rules at every level, you could just as well do:

result = ua[z, t] - uw[z, t] /. sol2
share|improve this answer
    
Thx for your answer –  Martin Wijaya Feb 5 '13 at 12:48
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.