0
$\begingroup$
EQU = {Y'[s] == -2 X[s] /Sin[ 2 Y[s]] Z[s], 
   X'[s] == Sin[X[s]] Cos[Y[s]]/3, Z'[s] == Cos[Y[s]] Cos[X[s]], 
   X[0] == 0, Y[0] == .3, Z[0] == 1, V1'[s] == Z[s]^2 Cos[Y[s]], 
   V1[0] == 0};
NDSolve[EQU, {X, Y, Z, V1}, {s, 0, 2}];
{x[u_], y[u_], z[u_], v1[u_]} = {X[u], Y[u], Z[u], V1[u]} /. First[%];
Plot[{x[s], y[s], z[s], v1[s]}, {s, 0, 2}, PlotLabel -> "Z X Y V1"]
v2[s_] = NIntegrate[z[s]^2 Cos[y[s]], {s, 0, 2}];
Plot[{x[s], y[s], z[s], v2[s]}, {s, 0, 2}, PlotLabel -> "Z X Y v2"]

EDIT1:

V1 is calculated inside NDSolve correctly but the same function (forgotten in the first instance) is calculated with a different name V2 outside NDSolve separately. I used earlier obtained function definitions for constituents as input in NIntegrate. But it does not work. There is something I miss in V2 ( where V1 calculated correctly from initial point of differential equation solution, but fails for V2.) Please help to find it out.

EDIT2

r1 = 3.25; ro = 1; lmax = 5;
EQN := {ph1'[
     l] == -Cos[ph1[l]]/
      rr1[l]*(2 rr1[l]^2 - 3 ro^2)/(rr1[l]^2 - ro^2), 
   rr1'[l] == Sin[ph1[l]], z1'[l] == Cos[ph1[l]], ph1[0] == 0, 
   zz1[0] == 0, rr1[0] == r1, vol1[l] == Pi rr1[l]^2 Cos[ph1[l]], 
   vol1[0] == 0};
NDSolve[EQN, {ph1, rr1, z1, vol1}, {l, 0, lmax}];
{r[u_], ph[u_], z[u_], 
   vol[u_]} := {rr1[u], ph1[u], z1[u], vol1[u]} /. First[%];
P1 = ParametricPlot[{z[l], r[l]}, {l, 0, lmax}, 
  GridLines -> Automatic, AspectRatio -> Automatic]
ParametricPlot[{z[l], ph[l]}, {l, 0, lmax}, AspectRatio -> 1]
volume[smax_] := NIntegrate[Pi r[l]^2 Cos[ph[l]], {l, 0, smax}]
Plot[{ph[l], volume[l]}, {l, 0, lmax}, AspectRatio -> 1]

The above does not work, although earlier it was ok.

EDIT3

Error messages before correct results suggested by MMM remedy for EDIT2:

enter image description here

$\endgroup$
2
  • 2
    $\begingroup$ Could you be MUCH clearer about what exactly your problem is? $\endgroup$
    – MarcoB
    Feb 7 '17 at 19:02
  • $\begingroup$ Hope it is more clear now. $\endgroup$
    – Narasimham
    Feb 7 '17 at 20:31
4
$\begingroup$

One way to do it is

{Y[s_], Z[s_], V1[s_]} = {Y[s], Z[s], V1[s]} /.First@NDSolve[EQU, {X, Y, Z, V1}, {s, 0, 2}];
V2[y_] := NIntegrate[Z[s]^2 Cos[Y[s]], {s, 0, y}];
Plot[{V1[s], V2[s]}, {s, 0, 2}, PlotRange -> All, 
 PlotStyle -> {Green, Directive[Red, Dashed]}, Frame -> True]

enter image description here

Please avoid using single capital alphabets.

Response to Edit 2

zz1[0]==0 should be z1[0]==0. vol1[l] should be vol1'[l].

r1 = 3.25; ro = 1; lmax = 5;
EQN = {ph1'[l] == -Cos[ph1[l]]/ rr1[l]*(2 rr1[l]^2 - 3 ro^2)/(rr1[l]^2 - ro^2), 
   rr1'[l] == Sin[ph1[l]], z1'[l] == Cos[ph1[l]], ph1[0] == 0, 
   z1[0] == 0, rr1[0] == r1, vol1'[l] == Pi rr1[l]^2 Cos[ph1[l]],vol1[0] == 0};
{rr1[l_], ph1[l_], z1[l_],vol1[l_]} = {rr1[l], ph1[l], z1[l], vol1[l]} /. 
   First[NDSolve[EQN, {ph1, rr1, z1, vol1}, {l, 0, lmax}]];
volme[x_] := NIntegrate[Pi rr1[l]^2 Cos[ph1[l]], {l, 0, x}]
Plot[{vol1[l], volme[l]}, {l, 0, 5}, AspectRatio -> 1]

enter image description here

$\endgroup$
3
  • $\begingroup$ It works but with error messages.. $\endgroup$
    – Narasimham
    Feb 9 '17 at 20:58
  • $\begingroup$ nice answer, note that you can extract the solution functions also like this: {X, Y, Z, V1} = NDSolveValue[EQU, {X, Y, Z, V1}, {s, 0, 2}]. $\endgroup$ Feb 9 '17 at 21:06
  • $\begingroup$ @AlbertRetey Thanks. I saw your two answers. I will try your idea. $\endgroup$
    – zhk
    Feb 10 '17 at 3:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.