2

In pstricks, psSolid and pstODEsolve the problem has been resolved using a function as "output-function" for a solution to an ode (refered to as the real ode) which is itself a solution to an ode, say an antiderivative. Many thanks!!! However I would like to compute not only one solution to the real ode but several say four or five. The following does indeed compute two data set but does not plot (even one) - I guess there might be more global variable -- but this is ways beyond my understanding

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\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 intermediateResult1.dat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\getFirstOfThree#1 #2 #3;{#1}
\def\getSecondOfThree#1 #2 #3;{#2}
\def\getThirdOfThree#1 #2 #3;{#3}

\IfFileExists{intermediateResult1.dat}{%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %  second pass:
  %
  %  post-process table --> z(second column)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % initialise Postscript variable to take the processed table
  \pstVerb{%
    true setglobal
    globaldict /processedTable1 {} put
    false setglobal
  }%
  \endlinechar=-1% suppress trailing space at input line end
  \newread\xIIfile%
  \openin\xIIfile=intermediateResult1.dat%
  %line-wise read table from `intermediateResult1.dat', treat second column
  \loop\read\xIIfile to \inputline%
  \ifeof\xIIfile\else%
   \pstODEsolve[algebraic,algebraicT]{processedVal1}{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 /processedTable1 [
         processedTable1
         \expandafter\getFirstOfThree\inputline;\space
         processedVal1 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 `intermediateResult1.dat'
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \pstODEsolve[algebraic,algebraicOutputFormat,saveData]{intermediateResult1}{%
      \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]))}%
}%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\IfFileExists{intermediateResult2.dat}{%
  \pstVerb{%
    true setglobal
    globaldict /processedTable2 {} put
    false setglobal
  }%
  \endlinechar=-1% suppress trailing space at input line end
  \newread\xIIIfile%
  \openin\xIIIfile=intermediateResult2.dat%
  \loop\read\xIIIfile to \inputline%
  \ifeof\xIIIfile\else%
   \pstODEsolve[algebraic,algebraicT]{processedVal2}{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)}%
     \pstVerb{
       true setglobal
       globaldict /processedTable2 [
         processedTable2
         \expandafter\getFirstOfThree\inputline;\space
         processedVal2 exch pop % insert second value and throw the first at t=0 away
         \expandafter\getThirdOfThree\inputline;] cvx put
       false setglobal
     }%
  \repeat%
}{%
    \pstODEsolve[algebraic,algebraicOutputFormat,saveData]{intermediateResult2}{%
      \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.7 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]))}%
}%


%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%

\begin{document}

%%
%%   The first solution should be plotted
%%
\IfFileExists{intermediateResult1.dat}{%
\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]{processedTable1}
  \axesIIID(1,4,1)(2,4.5,2)
\end{pspicture} 
}%
%%%%%%%%%%%%%%  
{
dummy text
}
\end{document}
  • So, if the particular problem of your previous question was solved, would you please click onto the check mark on the left side of the corresponding answer? (It is the only compensation we get for helping on this site.) – AlexG Jul 10 '15 at 8:31
  • Your ODE solution inside z(x) doesn't need 50 output points, since only the last one at the upper boundary is kept and plotted. So just take the mandatory minimum number of two (the integration interval borders) and throw away the unneeded one. I did this in the "second pass" branch. – AlexG Jul 10 '15 at 9:08
  • Some more comments: (1) All calculations can be done in the document preamble (before \begin{document}). (2) You must, of course, copy and paste the code block in the preamble as many times as you want different solutions of your main (real) ODE. Adjust the initial conditions in each copy. Replace processedTable in each code block with a unique name, as you did with intermediateResult. (3) Do you want to plot all computed trajectories into one figure? Then put all \listplotIIID into one pspicture environment. – AlexG Jul 10 '15 at 9:31
2

Again, as in the related question, we can apply the secondary ODE in z(x) onto the y coordinate of the surface vertices in a post-processing step.

Only one curve is plotted onto the surface. Adding more curves of the post-processed primary ODE solutions with modified initial conditions is left as an exercise. Just read the comments to your question and repeat the corresponding code blocks within sections "third pass" and "first pass". You will have to replace "intermediateResult" and "processedTable" in the repeated code with unique names, such as "intermediateResult1" and "processedTable1".

Now, in order to get the final picture, we need three passes of *.tex to *.pdf compilation. The second pass, which is processing the intermediate surface vertex coordinates, is quite greedy in terms of TeX memory. We use lualatex in DVI output mode here to avoid memory problems.

These are the three compilation passes for the code (example.tex) below to get the picture shown:

latex example
dvips example
ps2pdf -DNOSAFER example.ps

lualatex --output-format=dvi example
dvips example
ps2pdf -DNOSAFER example.ps

latex example
dvips example
ps2pdf -DNOSAFER example.ps

enter image description here


file example.tex:

\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 xyz.dat files (files with list of vertices)
\def\getFirstOfThree#1 #2 #3;{#1}
\def\getSecondOfThree#1 #2 #3;{#2}
\def\getThirdOfThree#1 #2 #3;{#3}

\IfFileExists{surf-sommets.dat}{%
  \endlinechar=-1% suppress trailing space at input line end
  \newread\xIIfile%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %  third pass
  %
  %  * post-process y coordinates in intermediateResult.dat (curve on surface)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  \pstVerb{% initialise Postscript variable to take the processed table
    true setglobal
    globaldict /processedTable {} put
    false setglobal
  }%
  \openin\xIIfile=intermediateResult.dat%
  %line-wise read table from `intermediateResult.dat', and treat second column: z(2nd 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 of pstODEsolve result and throw the first at t=0 away
        \expandafter\getThirdOfThree\inputline;] cvx put
      false setglobal
    }%
  \repeat%
}{
  \IfFileExists{intermediateSurf-sommets.dat}{%
    \endlinechar=-1% suppress trailing space at input line end
    \newread\xIIfile%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %  second pass
    %
    %  * post-process vertex y coordinates in intermediateSurf-sommets.dat
    %    write result into surf-sommets.dat
    %  * copy remaining intermediateSurf-*.dat to surf-*.dat files
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % post-process vertex y coordinates in intermediateSurf-sommets.dat
    \openin\xIIfile=intermediateSurf-sommets.dat%
    \pstVerb{ %ps2pdf will write surf-sommets.dat
      /surfsommets (surf-sommets.dat) (w) file def
      /numberstring 256 string def
      /num2str {numberstring cvs} def
    }
    %line-wise read table from `intermediateSurf-sommets.dat'
    \loop\read\xIIfile to \inputline%
    \ifeof\xIIfile
      \closein\xIIfile
      \pstVerb{surfsommets closefile}
    \else%
      %treat second column (output ODE)
      \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)}%
       %write resulting table into `surf-sommets.dat'
       \pstVerb{
         surfsommets \expandafter\getFirstOfThree\inputline; num2str writestring
         surfsommets ( ) writestring
         surfsommets processedVal exch pop num2str writestring
         surfsommets ( ) writestring
         surfsommets \expandafter\getThirdOfThree\inputline; num2str writestring
         surfsommets (\string\n) writestring
       }%
    \repeat%
    % copy remaining intermediateSurf-*.dat to surf-*.dat files
    \newwrite\copiedFile
    \openin\xIIfile=intermediateSurf-couleurs.dat%
    \immediate\openout\copiedFile=surf-couleurs.dat%
    \loop\read\xIIfile to \inputline%
    \ifeof\xIIfile
      \closein\xIIfile
      \closeout\copiedFile
    \else%
      \immediate\write\copiedFile{\inputline}
    \repeat%
    \openin\xIIfile=intermediateSurf-faces.dat%
    \immediate\openout\copiedFile=surf-faces.dat%
    \loop\read\xIIfile to \inputline%
    \ifeof\xIIfile
      \closein\xIIfile
      \closeout\copiedFile
    \else%
      \immediate\write\copiedFile{\inputline}
    \repeat%
    \openin\xIIfile=intermediateSurf-io.dat%
    \immediate\openout\copiedFile=surf-io.dat%
    \loop\read\xIIfile to \inputline%
    \ifeof\xIIfile
      \closein\xIIfile
      \closeout\copiedFile
    \else%
      \immediate\write\copiedFile{\inputline}
    \repeat%
  }{
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %  first pass
    %
    %  * solve main ODE (curve on surface) 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 in the second pass
      \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}

\begin{pspicture}(-0.7,-1.1)(3.2,1.15)%
  \psset[pst-solides3d]{viewpoint=20 5 20,Decran=20,lightsrc=viewpoint}%
  \IfFileExists{surf-sommets.dat}{%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %  third pass
    %
    %  plot readily processed surface from surf-*.dat files
    %  plot readily processed curve in PS variable `processedTable'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \psSolid[object=datfile, file=surf,
      base=0 3.14 0 6.28,
      inhue=0.01 0.1, hue=0.1 0.01,
      linewidth=0.5pt,
      grid=false,
      opacity=0.9,
      ngrid=200 70,
      intersectiontype=0,intersectionplan={[ 0 0 1 0] [0 1 0 0]},intersectioncolor=(blue) (red),
      intersectionlinewidth=1 2
    ]%
    \listplotIIID[linecolor=black,linewidth=1pt]{processedTable}
  }{%
    \IfFileExists{intermediateSurf-sommets.dat}{%
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %  second pass
      %
      %  do nothing here
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    }{%
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %  first pass
      %
      %  write the surface to intermediateSurf-*.dat files for postprocessing in the 2nd pass
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %the surface definition with the un-processed vertex coordinate y = t
      %(to be post-prcessed in the second pass)
      \defFunction[algebraic]{surf}(t,v){\r(t)*sin(v)}{t}{\r(t)*cos(v)}
      \psSolid[
        file=intermediateSurf,action=writesolid,
        object=surfaceparametree,
        function=surf,
        base=0 3.14 0 6.28,
        inhue=0.01 0.1, hue=0.1 0.01,
        linewidth=0.5pt,
        grid=false,
        opacity=0.9,
        ngrid=200 70,
        intersectiontype=0,intersectionplan={[ 0 0 1 0] [0 1 0 0]},intersectioncolor=(blue) (red),
        intersectionlinewidth=1 2
      ]%
    }%
  }%
  \axesIIID(1,4,1)(2,4.5,2)%
\end{pspicture}

\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.