Solve 1-D function that has 2-D interpolating functions

Question

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

Answer 1

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}]

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

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