# Using data from NDSolve into a secondary equation

I have this ODE,

f''''[y]+Z1*(6 f''[y]*f'''[y]*f'''[y]+3 f''[y]*f''[y]*f''''[y])- N1*N1*f''[y] == 0


subject to boundary conditions

f[h2] == F1/2, f[h1] == -(F1/2), f'[h1] == 0, f'[h2] == 0


where

h1 = -1 - m*x - a*Sin[2*\[Pi]*(x - t) + \[Phi]]
h2 = 1 + m*x + b*Sin[2*\[Pi]*(x - t)]

F1 = Q1 + a*Sin[2*\[Pi]*(x - t)]+b*Sin[2*\[Pi]*(x - t) + \[Phi]]
Z1 = 0.002; N1 = Sqrt[M1^2 + (1/k)]; m = 0.25; \[Phi] = 2 \[Pi]/3;k = 0.8;a = 0.3; b = 0.4;


The secondary equation is

SE = D[f[y], {y, 3}] + Z1*D[f[y], {y, 3}]*D[f[y], {y, 2}]*D[f[y], {y, 2}]+
2*D[f[y], {y, 2}]*D[f[y], {y, 2}]*D[f[y], {y, 2}]*D[f[y], {y, 3}]
- N1*N1*D[f[y], y]

1. I want to integrate the secondary equation over x=0..1 and t=0..1 and then plot it against Q1=-3..3 for three values of M1=1,2,3 with y -> 0.

2. How can I plot contours of f(x,y) for fix values of t=1 and M1=1?

• Zeta in Mathematica is the Riemann zeta function $\zeta(s)$. You mustn't use it as a constant. May 7, 2016 at 12:29
• @Artes thanks for the comment. I have made the changes.
– zhk
May 7, 2016 at 14:53
• Avoid using capital letters when assigning constants, Mathematica has an enormous library! May 7, 2016 at 15:17
• Where's the value of a and b? What's the value of y in SE? May 9, 2016 at 5:38
• @xzczd my bad. I have modified the post. Thanks!
– zhk
May 9, 2016 at 9:12

I am very scared of big equations. So whenever I see one, I go numeric (sometimes it is faster also).

m = 0.25; \[Phi] = 2 \[Pi]/3; k = 0.8; a = 0.3; b = 0.4;

h1 = -1 - m*x - a*Sin[2*\[Pi]*(x - t) + \[Phi]]
h2 = 1 + m*x + b*Sin[2*\[Pi]*(x - t)]

F1 = Q1 + a*Sin[2*\[Pi]*(x - t)] + b*Sin[2*\[Pi]*(x - t) + \[Phi]]
Z1 = 0.002; N1 = Sqrt[M1^2 + (1/k)];


I am going to show one example for

M1=1


You can use NDSolve for a given Q1, x, and t like

 Clear[f]
f[y_] = f[y] /.
Block[{Q1 = 1, x = 0.3, t = 0.2},
NDSolve[{f''''[y] + Z1*(6 f''[y]*f'''[y]*f'''[y] + 3 f''[y]*f''[y]*f''''[y]) -
N1*N1*f''[y] == 0, f[h2] == F1/2, f[h1] == -(F1/2), f'[h1] == 0, f'[h2] == 0},
f[y], y]];

SE = D[f[y], {y, 3}] + Z1*D[f[y], {y, 3}]*D[f[y], {y, 2}]*D[f[y], {y, 2}] +
2*D[f[y], {y, 2}]*D[f[y], {y, 2}]*D[f[y], {y, 2}]*D[f[y], {y, 3}] -
N1*N1*D[f[y], y];

SE /. y -> 0


-2.44329

You can create a {x,t,SE} and use any numerical method to integrate. A lazy person like me will simply use Sum (which works nice for small dx and dt)

data= Table[
{Q,

Sum[
Clear[f];
f[y_] =
f[y] /. Block[{Q1 = Q, x = x1, t = t1},
NDSolve[{f''''[y] +
Z1*(6 f''[y]*f'''[y]*f'''[y] + 3 f''[y]*f''[y]*f''''[y]) -
N1*N1*f''[y] == 0,
f[h2] == F1/2, f[h1] == -(F1/2), f'[h1] == 0, f'[h2] == 0},
f[y], y]][[1]];

SE = D[f[y], {y, 3}] +
Z1*D[f[y], {y, 3}]*D[f[y], {y, 2}]*D[f[y], {y, 2}] +
2*D[f[y], {y, 2}]*D[f[y], {y, 2}]*D[f[y], {y, 2}]*
D[f[y], {y, 3}] - N1*N1*D[f[y], y];

SE /. y -> 0,

{x1, 0., 1, 0.1}, {t1, 0., 1, 0.1}]}

, {Q, -3., 3., 0.5}]


{{-3., 1134.48}, {-2.5, 855.009}, {-2., 639.126}, {-1.5, 461.605}, {-1., 305.008}, {-0.5, 157.953}, {0., 13.8868}, {0.5, -129.68}, {1., -271.471}, {1.5, -406.275}, {2., \ -524.337}, {2.5, -610.275}, {3., -641.51}}

And then a aimple ListlinePlot

ListLinePlot[data, Mesh -> All]


For the ContourPlot, you can fix a value for Q1 and generate a Table in a same way.

M1 = 1; Q0 = 1; t0 = 1;
data = Table[
Clear[f];
f[y_] =
f[y] /. Block[{Q1 = Q0, x = x1, t = t0},
NDSolve[{f''''[y] +
Z1*(6 f''[y]*f'''[y]*f'''[y] + 3 f''[y]*f''[y]*f''''[y]) -
N1*N1*f''[y] == 0,
f[h2] == F1/2, f[h1] == -(F1/2), f'[h1] == 0, f'[h2] == 0},
f[y], y]][[1]];

{x1, y1, f[y1]},

{x1, 0., 1, 0.1}, {y1, 0., 1, 0.1}];
data = Flatten[data, 1];

ListContourPlot[data]


For a better result and smoother plot you can reduce the stepsize. I use 0.1. You can go further below like 0.01 or less (I am too impatient for that). Do some trial value and you can see if your result is varying much due to the change and you can guess the optimum step size.

For multiple variable

a = 0.3; b = 0.4; m = 0.25; \[Phi] = 2 \[Pi]/3; M1 = 1; k = 0.8;
Nt = 0.8;Nb = 0.4; Pr = 0.2; L = 0.1; Bm = 4; Bh = 2; ;Gr = 0.7;
Br = 1; \[Zeta] = 0.002; Rn = 0.6;

Clear[f1, f2, f3]
{f1[y_], f2[y_], f3[y_]} = {f1[y], f2[y], f3[y]} /.
Block[{Q1 = 1, x = 0.3, t = 0.2, N1 = 0.1, F1 = 0.3},
NDSolve[{f1''''[
y] + \[Zeta]*(6 f1''[y]*f1'''[y]*f1'''[y] +
3 f1''[y]*f1''[y]*f1''''[y]) - N1*N1*f1''[y] + Gr*f2'[y] +
Br*f3'[y] == 0,
(1 + 1)*f2''[y] + f3'[y]*f2'[y] + f2'[y]*f2'[y] == 0,
f3''[y] + f2''[y] == 0,
f1[h2] == F1/2, f1[h1] == -(F1/2), f1'[h1] == 0,
f1'[h2] == 0, f3'[h1] == f3[h1], f3'[h2] == (1 - f3[h2]),
f2'[h1] == Bh*f2[h1], f2'[h2] == (1 - f2[h2])}, {f1[y], f2[y],
f3[y]},
y]][[1]];

{f1[0], f2[0], f3[0]}


{-0.0191435, 0.525077, 0.379438}

• Thanks dear for your efforts. What if we have more than one dependent variables? How we can then write f[y_] =f[y] /. for it? It seems that the rest of the structure will remain the same.
– zhk
May 9, 2016 at 12:34
• Just like that. f[x1_,x2_,x3_...] = f[x1,x2,x3...]/.. Make sure the same functional form apears at the end of your NDSolve. I would suggest using a region for the NDSolve. I mean instead of using f[y], y], use f[y], {y,-1.,1.}]. Since the result comes in form of interpolation, this might give you a better (and faster) result. May 9, 2016 at 14:57
• I think you miss understood me. If I have three dependent variables say f1[y], f2[y], f3[y] then how we will write the functional for this? This is the case when the secondary equation contains all of the three f1,f2,f3?
– zhk
May 9, 2016 at 16:23
• Got it. Again, just like that. {f1[y_], f2[y_], f3[y_]} = {f1[y], f2[y], f3[y]}/.. May 9, 2016 at 16:59
• I am getting this error  Lists {f1[y_],f2[y_],f3[y_]} and {{<<1>>,<<1>>,<<1>>}} are not the same shape.
– zhk
May 10, 2016 at 2:20