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 need some helps on NDSolve. Here are my equations:

k1 = 0.000192;
k2 = 0.000159;

Clear[x, y, z, h]
eqnf1 = {x'[t] == -k1*x[t]*y[t]/h[t] + k2*z[t],
   y'[t] == 3*(-k1*x[t]*y[t]/h[t] + k2*z[t]),
   z'[t] == k1*x[t]*y[t]/h[t] - k2*z[t],
   h'[t] == -3*(-k1*x[t]*y[t]/h[t] + k2*z[t]),
   x[0] == .001, z[0] == 0, y[0] == 0.2, h[0] == .001};

sol2 = NDSolve[eqnf1, {x, y, z, h}, {t, 0, 1000}][[1]]

X[t_] = x[t] /. sol2;
Y[t_] = y[t] /. sol2;
Z[t_] = z[t] /. sol2;
H[t_] = h[t] /. sol2;

I used this part to solve y[t], which will be using later for my other function. And here is 2nd part of my function:

Da = 1.33*10 - 4;
eqnm1 = {D[yc[x, t], {t, 1}] == 
    Da* D[ D[yc[x, t], {x, 1}], {x, 1}] ,
   yc[x, 0] == 0, 
   yc[0, t] == Z[t], (D[yc[x, t], {x, 1}] /. x -> 1) == 0};
solmem1 = 
 First[yc /. NDSolve[eqnm1, yc, {x, 10^-20, 1}, {t, 0, 1000}]]

Here I obtained another interpolating function, yc[x,t], which will be used in the last part of the problem that I'm not able to solve. The last part of the problem needs function from the 1st and 2nd parts - y[t] & yc[x,t] which both are interpolating functions. And the 3rd part is only 1D function, which is xs[t]. I tried

eqns1 = {xs'[t] == -k1*xs[t]*First[Evaluate[y[t] /. eqnf1]]/hs[t] + 
     k2*First[Evaluate[yc[x, t] /. x -> 1 /. solmem1]],
   hs'[t] == -3*(-k1*xs[t]*First[Evaluate[y[t] /. eqnf1]]/hs[t] + 
       k2*First[Evaluate[yc[x, t] /. x -> 1 /. solmem1]]),
   xs[0] == 0, hs[0] == 1};

But it has some errors that I think is caused by the 2nd part of the problem that is 2-D (ys[x,t]). Does anyone know how to solve xs(t) and hs(s)? TThanks

share|improve this question
    
I suspect you want ... /. sol2 (rather than eqnf1) and also ... /. ss1 needs ss1 defined. Probably that was to come from the modification solmem1 = First[yc /. (ss1 = NDSolve[eqnm1, yc, {x, 10^-20, 1}, {t, 0, 1000}])] –  Daniel Lichtblau Dec 13 '12 at 0:32
    
Yeah sorry... I had some typoe. Hopefully it's all corrected now –  DumbleKo Dec 13 '12 at 0:43
    
No, not fixed. Your replacements are not going to be of much use. (1) What are you expecting to obtain when you do ... /. eqnf1 ? (Hint: Per my earlier comment, you possibly want /. sol2 in such places.) (2) What do you expect to obtain when you do ... /. solmem1? (Hint: look at how I defined ss1 in my earlier comment. Second hint: By "look at..." I mean specifically "Do not ignore...") All this might not get you what you want but at least it will get you using viable replacement rules. Also it gives something that I can confirm NDSolve will handle. –  Daniel Lichtblau Dec 13 '12 at 14:52
    
Hi Daniel, Thanks for you help. Next time I'll check my equations more carefully. –  DumbleKo Dec 13 '12 at 23:09

1 Answer 1

up vote 3 down vote accepted

I'll redo this from the top, simplifying a few things in your systems. Also I replace capitalized variables by ones starting in lower case; that's good Mathematica practice.

k1 = 0.000192;
k2 = 0.000159;

eqnf1 = {x'[t] == -k1*x[t]*y[t]/h[t] + k2*z[t], y'[t] == 3*x'[t], 
  z'[t] == -x'[t], h'[t] == -3*x'[t], x[0] == .001, z[0] == 0, 
  y[0] == 0.2, h[0] == .001}; sol2 = 
 NDSolve[eqnf1, {x, y, z, h}, {t, 0, 1000}][[1]];

xX[t_] = x[t] /. sol2;
yY[t_] = y[t] /. sol2;
zZ[t_] = z[t] /. sol2; (* We only use this one later *)
hH[t_] = h[t] /. sol2;

Da = 1.33*10 - 4;
eqnm1 = {D[yc[x, t], {t, 1}] == 
    Da*D[D[yc[x, t], {x, 1}], {x, 1}] + k2*yc[x, t] - 1*yc[x, t], 
   yc[x, 0] == 0, 
   yc[0, t] == zZ[t], (D[yc[x, t], {x, 1}] /. x -> 1) == 0};
solmem1 = 
 First[yc /. (ss1 = NDSolve[eqnm1, yc, {x, 10^-20, 1}, {t, 0, 1000}])]

Here are corrected replacements for eqns1. Again I simplify by avoiding repetition of a long rhs.

eqns1 = {xs'[t] == -k1*xs[t]*First[Evaluate[y[t] /. sol2]]/hs[t] + 
     k2*First[Evaluate[yc[x, t] /. x -> 1 /. ss1]], 
   hs'[t] == -3*xs'[t], xs[0] == 0, hs[0] == 1};

We can solve for this just fine.

ss2 = NDSolve[eqns1, {xs, hs}, {t, 0, 1000}][[1]]

Not sure what you are expecting. But here is what you get.

Plot[{xs[t]} /. ss2, {t, 0, 1000}]

enter image description here

Plot[hs[t] /. ss2, {t, 0, 1000}]

enter image description here

Final remark: Next time please read the comment more carefully. That would save us both time and work.

share|improve this answer

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.