2

When using pstODEsolve with the option algebraicOutputFormat I would like to use a function \z which itself is defined by the "final value" of a solution to another ode and the argument #1 being the time to be integrated over

\documentclass[border=0pt,10pt,dvipsnames]{standalone}
\usepackage{amsmath}
\usepackage{pstricks}
\usepackage{dashrule}
\usepackage{pst-plot}%, pst-math}
\usepackage{pst-func}
\usepackage{pst-ode}
\usepackage{pst-3dplot}
\usepackage{pst-solides3d}
\usepackage[off]{auto-pst-pdf}
\usepackage[nomessages]{fp}
\psset{algebraic,unit=2}
\usepackage[active,tightpage]{preview}
\PreviewBorder=12pt
\PreviewEnvironment{pspicture}


\makeatletter
\define@key[psset]{}{valuerange}[1-]{\expandafter\pst@getrange#1\@nil}
\def\pst@getrange#1-#2\@nil{%
  \ifx\relax#1\relax \def\pst@startvalue{1 }\else\def\pst@startvalue{#1 }\fi%
  \ifx\relax#2\relax \def\pst@endvalue{1e32 }\else\def\pst@endvalue{#2 }\fi%
}
\psset{valuerange=1-}

\def\listplotIIID{\def\pst@par{}\pst@object{listplotIIID}}
\def\listplotIIID@i#1{%
\@nameuse{beginplot@\psplotstyle}%
\addto@pscode{%
   /viewpointXYZ {\pst@solides@viewpoint} def
   /Decran \pst@solides@Decran\space def % distance de l'ecran
    viewpointXYZ /ZpointVue ED /YpointVue ED /XpointVue ED
  /THETA {YpointVue XpointVue atan} bind def
  /PHI   {ZpointVue XpointVue dup mul YpointVue dup mul add sqrt atan} bind def
  /Dobs  {XpointVue dup mul YpointVue dup mul add ZpointVue dup mul add sqrt} bind def
  XpointVue YpointVue ZpointVue /viewpoint defpoint3d
    /XYZ [#1] def
  /@tabXYZ [
    0 3 XYZ length 3 sub {
     /i exch def
     XYZ i get
     XYZ i 1 add get
     XYZ i 2 add get
     i 3 div dup /ii ED  
       \pst@startvalue ge {
         ii \pst@endvalue le { 
         3dto2d
         \pst@number\psunit mul exch
         \pst@number\psunit mul exch
       } { pop pop pop } ifelse
     }{ pop pop pop } ifelse
    }  for
    ] bind def
 [ @tabXYZ aload pop
    }%
\@nameuse{endplot@\psplotstyle}}%
\makeatother



\def\r(#1){sin(#1)}
\def\rs(#1){cos(#1)}
\def\rss(#1){-sin(#1)}
\def\rssrs(#1){-sin(#1)*cos(#1)}
\def\zsq(#1){sin(#1)*sin(#1)*(1+cos(#1)+sin(2*#1)*sin(2*#1)/16)}
\def\zsszs(#1){(1/2)*sin(2*#1)*(1+cos(#1)+sin(2*#1)*sin(2*#1)/16)+(1/2)*sin(#1)*sin(#1)*(-sin(#1)+sin(4*#1)/8)}
%%%%%
%%%%%   How to define a function as the "last" value of a numerical solution
%%%%%   to an ODE where the variable #1 is the time to be integrated over
%%%%%   starting from time = 0 in this example
%%%%%
%%%%%   the following is of course not defined!
%%%%%

\def\z(#1){ "last value of" \pstODEsolve[algebraic,algebraicT]{test}{1}{0}{#1}{50}{0 1}{ 1 | 
 -sin(x[0])*sqrt(1+cos(x[0])+sin(2*x[0])*sin(2*x[0])/16)} }


\begin{document}
\begin{pspicture}(-1.2,-1)(3.1,1.1)
\psset[pst-solides3d]{viewpoint=20 5 20,Decran=20,lightsrc=viewpoint}
\pstODEsolve[algebraic,algebraicOutputFormat]{test}{\r(x[2])*sin(x[0]) | \z(x[2]) 
| \r(x[2])*cos(x[0])}{0}{22}{200}{1.57 1 0.3 0}{ x[1] | -2*\rs(x[2])*x[1]*x[3]/(\r(x[2]))| x[3] | -(\rssrs(x[2])+\zsszs(x[2]))*x[3]*x[3]/(\rs(x[2])*\rs(x[2])+\zsq(x[2]))
+\rs(x[2])*\r(x[2])*x[1]*x[1]/(\rs(x[2])*\rs(x[2])+\zsq(x[2]))}
\listplotIIID[linecolor=black,linewidth=1pt]{test}
\axesIIID(1,4,1)(2,4.5,2)
\end{pspicture}
\end{document}
  • I am afraid you cannot nest \pstODEsolve in that way. – AlexG Jul 8 '15 at 10:56
  • I am going to try something, please be patient. – AlexG Jul 8 '15 at 12:15
  • Something is ready, hope it is what you want. – AlexG Jul 8 '15 at 14:46
3

\pstODEsolve cannot be nested.

Fortunately, in the present case, we can work-around the problem, since the requested function z(x) is to be applied as a post-processing step on the second column of the solution of your main ODE.

This can be achieved in a two-step procedure.

1) solution of the main ODE without applying z(x) on the second element of the output vector:

\pstODEsolve[algebraic,algebraicOutputFormat,saveData]{intermediateResult}{%
  \r(x[2])*sin(x[0]) |
  x[2] | % we process this column later
  \r(x[2])*cos(x[0])
}{0}{22}{200}{1.57 1 0.3 0}{
 x[1] | -2*\rs(x[2])*x[1]*x[3]/(\r(x[2]))| x[3] |
 -(\rssrs(x[2])+\zsszs(x[2]))*x[3]*x[3]/(\rs(x[2])*\rs(x[2])+\zsq(x[2]))
 +\rs(x[2])*\r(x[2])*x[1]*x[1]/(\rs(x[2])*\rs(x[2])+\zsq(x[2]))
}

and write the solution table into a file intermediateResult.dat

2) Read the table from the intermediate file in a line-wise fashion, apply the other ODE

                                        % two output points--v
\pstODEsolve[algebraic,algebraicT]{processedValue}{1}{0}{#1}{2}{0 1}{
  1 | -sin(x[0])*sqrt(1+cos(x[0])+sin(2*x[0])*sin(2*x[0])/16)
}%

on the second element of the input line, take the last value from the result stored in processedValue and reassemble the table.


Since the PS --> PDF conversion step involves writing and reading text files, ps2pdf must be given the option -dNOSAFER.

Run: latex <...> dvips <...> ps2pdf -dNOSAFER <...>.ps

twice on the following input file:

\documentclass[border=10pt]{standalone}
\usepackage{pstricks}
\usepackage{pst-plot}%, pst-math}
\usepackage{pst-ode}
\usepackage{pst-solides3d}
\psset{algebraic,unit=2}

\makeatletter
\define@key[psset]{}{valuerange}[1-]{\expandafter\pst@getrange#1\@nil}
\def\pst@getrange#1-#2\@nil{%
  \ifx\relax#1\relax \def\pst@startvalue{1 }\else\def\pst@startvalue{#1 }\fi%
  \ifx\relax#2\relax \def\pst@endvalue{1e32 }\else\def\pst@endvalue{#2 }\fi%
}
\psset{valuerange=1-}

\def\listplotIIID{\def\pst@par{}\pst@object{listplotIIID}}
\def\listplotIIID@i#1{%
\@nameuse{beginplot@\psplotstyle}%
\addto@pscode{%
   /viewpointXYZ {\pst@solides@viewpoint} def
   /Decran \pst@solides@Decran\space def % distance de l'ecran
    viewpointXYZ /ZpointVue ED /YpointVue ED /XpointVue ED
  /THETA {YpointVue XpointVue atan} bind def
  /PHI   {ZpointVue XpointVue dup mul YpointVue dup mul add sqrt atan} bind def
  /Dobs  {XpointVue dup mul YpointVue dup mul add ZpointVue dup mul add sqrt} bind def
  XpointVue YpointVue ZpointVue /viewpoint defpoint3d
    /XYZ [#1] def
  /@tabXYZ [
    0 3 XYZ length 3 sub {
     /i exch def
     XYZ i get
     XYZ i 1 add get
     XYZ i 2 add get
     i 3 div dup /ii ED  
       \pst@startvalue ge {
         ii \pst@endvalue le { 
         3dto2d
         \pst@number\psunit mul exch
         \pst@number\psunit mul exch
       } { pop pop pop } ifelse
     }{ pop pop pop } ifelse
    }  for
    ] bind def
 [ @tabXYZ aload pop
    }%
\@nameuse{endplot@\psplotstyle}}%
\makeatother

\def\r(#1){sin(#1)}
\def\rs(#1){cos(#1)}
\def\rss(#1){-sin(#1)}
\def\rssrs(#1){-sin(#1)*cos(#1)}
\def\zsq(#1){sin(#1)*sin(#1)*(1+cos(#1)+sin(2*#1)*sin(2*#1)/16)}
\def\zsszs(#1){(1/2)*sin(2*#1)*(1+cos(#1)+sin(2*#1)*sin(2*#1)/16)+(1/2)*sin(#1)*sin(#1)*(-sin(#1)+sin(4*#1)/8)}

%for parsing lines of intermediateResult.dat
\def\getFirstOfThree#1 #2 #3;{#1}
\def\getSecondOfThree#1 #2 #3;{#2}
\def\getThirdOfThree#1 #2 #3;{#3}

\IfFileExists{intermediateResult.dat}{%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %  second pass:
  %
  %  post-process table --> z(second column)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % initialise Postscript variable to take the processed table
  \pstVerb{%
    true setglobal
    globaldict /processedTable {} put
    false setglobal
  }%
  \endlinechar=-1% suppress trailing space at input line end
  \newread\xIIfile%
  \openin\xIIfile=intermediateResult.dat%
  %line-wise read table from `intermediateResult.dat', treat second column
  \loop\read\xIIfile to \inputline%
  \ifeof\xIIfile\else%
    \pstODEsolve[algebraic,algebraicT]{processedVal}{1}{0}{\expandafter\getSecondOfThree\inputline;}{2}{0 1}{ 1 |
     -sin(x[0])*sqrt(1+cos(x[0])+sin(2*x[0])*sin(2*x[0])/16)}%
     %process input line and append result to postscript variable `processedTable' (executable list)
     \pstVerb{
       true setglobal
       globaldict /processedTable [
         processedTable
         \expandafter\getFirstOfThree\inputline;\space
         processedVal exch pop % insert second value and throw the first at t=0 away
         \expandafter\getThirdOfThree\inputline;] cvx put
       false setglobal
     }%
  \repeat%
}{%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %  first pass
  %
  %  solve main ODE without special treatment for second column of the solution table
  %  write solution to file `intermediateResult.dat'
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \pstODEsolve[algebraic,algebraicOutputFormat,saveData]{intermediateResult}{%
      \r(x[2])*sin(x[0]) |
      x[2] | % we process this column later
      \r(x[2])*cos(x[0])
    }{0}{22}{200}{1.57 1 0.3 0}{
     x[1] | -2*\rs(x[2])*x[1]*x[3]/(\r(x[2]))| x[3] |
     -(\rssrs(x[2])+\zsszs(x[2]))*x[3]*x[3]/(\rs(x[2])*\rs(x[2])+\zsq(x[2]))
     +\rs(x[2])*\r(x[2])*x[1]*x[1]/(\rs(x[2])*\rs(x[2])+\zsq(x[2]))}%
}%

\begin{document}

\IfFileExists{intermediateResult.dat}{%
  %plot the resulting graph
  \begin{pspicture}(-1.2,-1)(3.1,1.1)
  \psset[pst-solides3d]{viewpoint=20 5 20,Decran=20,lightsrc=viewpoint}
  \listplotIIID[linecolor=black,linewidth=1pt]{processedTable}
  \axesIIID(1,4,1)(2,4.5,2)
  \end{pspicture}
}{
  dummy text
}  

\end{document}

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.