How manipulate NDSolve Interpolating Functions to be plotted or reused in further symbolic expressions

Question

I've been working with these non-linear differential equations. I was trying to get a symbolic solution using DSolve, however I only succeed in solving the system using NDSolve. The question here is, how can I reuse the obtained NDSolve interpolating functions, Y and Z, to obtain the derivatives of the expressions: FelastY and FelastZ (shown below) and then plot them to see their behavior as functions of theta. How can I manipulate or reuse numerical solutions to plot symbolic expressions? I only get an error message that says: "ReplaceAll::reps:{sol} is neither a list of replacement rules nor a valid dispatch table." Hope you can help me.

(* Parameters: *)

β = Pi;
c = 2;
m = 5;
exc = 26 10 - 4;
Ω = 1;
g = 9.81;
d = 0.0254;
ad = 0.3;
a = ad*d;
young = 210 10^9;
iner = (Pi/64)*d^4;
ld = 24;
L = ld*d;
kξ = 6.5;
kη = 7.5;
kξη = kξ/5;
ψ = ArcCos[(Y[θ]*Cos[θ] + 
  Z[θ]*Sin[θ])/(Sqrt[
  Y[θ]^2 + Z[θ]^2])];
G = (1 + Sin[ψ])/2;
c2t = TrigExpand[Cos[2*θ]];
s2t = TrigExpand[Sin[2*θ]];
ctb = TrigExpand[Cos[θ - β]];
stb = TrigExpand[Sin[θ - β]];

(* Equations: *)

eq1 = (m*Y''[θ]) + (c*Y'[θ]) + (kξ + kη)*
 Y[θ]/2 + 
G/2*(((kξ - kη)*c2t) - (2*kξη*s2t))*
 Y[θ] + 
G/2*(((kξ - kη)*s2t) + (2*kξη*c2t))*
 Z[θ] == (m*exc*(Ω^2)*ctb) - (m*g);

eq2 = (m*Z''[θ]) + (c*Z'[θ]) + (kξ + kη)*
 Z[θ]/2 + 
G/2*(((kξ - kη)*c2t) - (2*kξη*s2t))*
 Z[θ] + 
G/2*(((kξ - kη)*s2t) + (2*kξη*c2t))*
 Y[θ] == (m*exc*(Ω^2)*stb);

sol2 = NDSolve[{eq1, eq2, Y[0] == -7, Z[0] == 0.1, Y'[0] == 0, 
Z'[0] == 0}, {Y[θ], Z[θ]}, {θ, 0, 200*Pi}]

g4 = ParametricPlot[
Evaluate[{Y[θ], Z[θ]} /. sol2], {θ, 0, 200*Pi},
PlotRange -> All]

These are the expressions in which I would like to reuse the numeric solutions Y and Z (interpolating functions) that I obtained from NDSolve:

FelastY[Y_?NumericQ, Z_?NumericQ] := (kξ + kη)*Y[θ]/2 + 
G/2*(((kξ - kη)*c2t) - (2*kξη*s2t))*Y[θ] + 
G/2*(((kξ - kη)*s2t) + (2*kξη*c2t))*Z[θ]

FelastY[Y[th] /. sol2]

FelastZ[Y_?NumericQ, Z_?NumericQ] := (kξ + kη)*Z[θ]/2 + 
G/2*(((kξ - kη)*c2t) - (2*kξη*s2t))*Z[θ] + 
G/2*(((kξ - kη)*s2t) + (2*kξη*c2t))*Y[θ]

FelastZ[Z[th] /. sol2]

Answer 1

Your question is a combination of simple mistakes, including:

1. Confounding function relationship (in your case, Y and Z) with function value (in your case, Y[th] and Z[th]). Remember these 2 concepts are strictly identified in Mathematica.

2. Abusing ?NumericQ. I guess you add it to your code simply because you found some questions about differential equations are solved with it, but, please don't use it if you don't understand its meaning, actually it's not related to your problem at all.

3. Abusing SetDelayed. For this part, just think about why x1 = x2; f[x2_] := 2 x1; f[3] doesn't output 6.

4. sol2 is blue i.e. the numeric solution isn't even assigned to sol2.

Anyway here is the fixed code.

Code lines before NDSolve[……] are unchanged so they're omitted in this answer.

{solY, solZ} = 
  NDSolveValue[{eq1, eq2, Y[0] == -7, Z[0] == 0.1, Y'[0] == 0, Z'[0] == 0}, {Y, 
    Z}, {θ, 0, 200*Pi}, MaxSteps -> Infinity];


g4 = ParametricPlot[{solY[θ], solZ[θ]}, {θ, 0, 200*Pi}, 
  PlotRange -> All, PlotPoints -> 100]

FelastY[Y_, Z_] = (kξ + kη)*Y[θ]/2 + 
   G/2*(((kξ - kη)*c2t) - (2*kξη*s2t))*Y[θ] + 
   G/2*(((kξ - kη)*s2t) + (2*kξη*c2t))*Z[θ];


FelastZ[Y_, Z_] = (kξ + kη)*Z[θ]/2 + 
   G/2*(((kξ - kη)*c2t) - (2*kξη*s2t))*Z[θ] + 
   G/2*(((kξ - kη)*s2t) + (2*kξη*c2t))*Y[θ];

ParametricPlot[{FelastY[solY, solZ], FelastZ[solY, solZ]}, {θ, 0, 200*Pi}, 
 PlotPoints -> 100]