home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 509.lha / ATS_v2.0 / ats.cp next >
Text File  |  1991-05-06  |  12KB  |  669 lines
  1. //        ret.terms[i]=complex(0.0,0.0);
  2. //
  3. //    C++ Routines to manipulate Taylor series for solving a differential
  4. //  equation.
  5. //
  6. //  By Dennis J. Darland 
  7. //  4/3/91
  8. //  Public Domain
  9. //  Version 1.6
  10. //  Many enhancements expected.
  11. //  PROVIDED WITH NO WARRANTY.
  12. //  
  13. #include "ats.h"
  14. t_s::t_s()
  15.     {
  16.     alloc_terms=MAX_TERMS;
  17.     }
  18. t_s::~t_s()
  19.     {
  20.     }
  21. t_s::t_s(t_s &a)
  22.     {
  23.     alloc_terms=MAX_TERMS;
  24.     set_first(a.first_term);
  25.     set_last(a.last_term);
  26.     set_x0(a.x0);
  27.     set_h(a.h);
  28.     set_terms(*this,a.terms);
  29.     }
  30. void t_s::set_first(int i)
  31.     {
  32.     first_term=i;
  33.     }
  34. void t_s::set_last(int i)
  35.     {
  36.     last_term=i;
  37.     }
  38. void t_s::set_h(complex h_in)
  39.     {
  40.     h=h_in;
  41.     }
  42. void t_s::set_x0(complex x0_in)
  43.     {
  44.     x0=x0_in;
  45.     }
  46. void set_data(t_s &a,t_s &b)
  47.     {
  48.     a.set_first(b.first_term);
  49.     a.set_last(b.last_term);
  50.     a.set_x0(b.x0);
  51.     a.set_h(b.h);
  52.     a.alloc_terms=MAX_TERMS;
  53.     }
  54. void set_terms(t_s &ret,complex *values)
  55.     {
  56.     register int i;
  57.     
  58.     for (i=0;i<MAX_TERMS;i++)
  59.         ret.terms[i]=values[i];
  60.     }
  61. complex ats(int m,t_s &a,t_s &b,int j)
  62.     {
  63.     // function to implement Liebniz rule
  64.     complex cret;
  65.     register int i;
  66.     
  67.     cret=complex(0.0,0.0);
  68.     
  69.     for (i=j;i<=m;i++)
  70.         {
  71.         cret += a.terms[i]*b.terms[m-i];
  72.         }
  73.     return(cret);
  74.     }
  75.  
  76. complex att(int m,t_s &a,t_s &b,int j)
  77.     {
  78.     // function for modified Liebniz rule
  79.     complex ret;
  80.     register int i;
  81.     
  82.     ret=complex(0.0,0.0);
  83.     if (j > m)
  84.         return(ret);    
  85.     for (i=j;i<=m;i++)
  86.         {
  87.         ret += a.terms[i]*b.terms[m-i+1]*complex(double(m-i+1),0.0);
  88.         }
  89.     ret /= complex(double(m+1),0.0);
  90.     return(ret);
  91.     }
  92. t_s & t_s::operator+=(t_s &b)
  93.     {
  94.     // Add taylor series
  95.     register int i;
  96.     
  97.     if ((first_term != b.first_term)
  98.     || (last_term != b.last_term))
  99.         report_error(0);
  100.     set_data(*this,b);
  101.     for (i=first_term; i < last_term;i++)
  102.         terms[i]+=b.terms[i];
  103.     return (*this);
  104.     }
  105. t_s & t_s::operator-=(t_s &b)
  106.     {
  107.     // Subtract taylor series
  108.     register int i;
  109.     
  110.     if ((first_term != b.first_term)
  111.     || (last_term != b.last_term))
  112.         report_error(0);
  113.     set_data(*this,b);
  114.     for (i=first_term; i < last_term;i++)
  115.         terms[i]-=b.terms[i];
  116.     return (*this);
  117.     }
  118. t_s& t_s::operator=(t_s &b)
  119.     {
  120.     // assign taylor series
  121.     register int i;
  122.     set_data(*this,b);
  123.     for (i=0; i < b.last_term;i++)
  124.         this->terms[i]=b.terms[i];
  125.     return(*this);
  126.     }
  127. void display(t_s &a,int plot)
  128.     {
  129.     static pt=0;
  130.     FILE *plot_cmds;
  131.     char fname[64];
  132.     char cmd[81];
  133.     
  134.     // display taylor series
  135.     register int i;
  136.     
  137.     sprintf(fname,"ram:mapleplot_%d",pt);
  138.     pt++;
  139.     plot_cmds = fopen("ram:plottmp","w");
  140.     fprintf(plot_cmds,"plotdevice := amiga;\n");
  141.     fprintf(plot_cmds,"plot([");
  142.     cout << "alloc_terms=" << a.alloc_terms << "\n"; 
  143.     cout << "first_term =" << a.first_term << "\n"; 
  144.     cout << "last_term=" << a.last_term << "\n";
  145.     cout << "h=" << a.h << "\n";
  146.     cout << "x0=" << a.x0 << "\n";
  147.     for (i=0;i<a.last_term;i++)
  148.         {
  149.         cout << "terms[" << i << "]= " << real(a.terms[i]) << " + "
  150.         << imag(a.terms[i]) << " I\n";
  151.         if (i != 0)
  152.             fprintf(plot_cmds,",");
  153.         fprintf(plot_cmds,"%d,sign(%g)*(ln(1+abs(%g)))\n",
  154.         i,real(a.terms[i]),real(a.terms[i]));
  155.         }
  156.     fprintf(plot_cmds,"]);\n");
  157.     fprintf(plot_cmds,"quit\n");
  158.     fclose(plot_cmds);
  159.     sprintf(cmd,"iconz <ram:plottmp >%s pltcvt",fname);
  160.     system(cmd);
  161.     if (plot>1)
  162.         {
  163.         sprintf(cmd,"maple <%s",fname);
  164.         system(cmd);
  165.         }
  166.     }
  167. void t_s::report_error(int no)
  168.     {
  169.     cout << "ats error = " << no;
  170.     }
  171. t_s operator+(t_s &a,t_s &b)
  172.     {
  173.     // Add two taylor series
  174.     register int i;
  175.     t_s ret;
  176.     complex ac,bc;
  177.     
  178.     set_data(ret,a);
  179.     for (i=a.first_term;i<a.last_term;i++)
  180.         {
  181.         ac=a.terms[i];
  182.         bc=b.terms[i];
  183.         ret.terms[i]=ac+bc;
  184.         }
  185.     return(ret);
  186.     }
  187.  
  188. t_s operator-(t_s &a,t_s &b)
  189.     {
  190.     // Subtract two taylor series
  191.     register int i;
  192.     t_s ret;
  193.     complex ac,bc;
  194.         
  195.     set_data(ret,a);
  196.     
  197.     for (i=a.first_term;i<a.last_term;i++)
  198.         {
  199.         ac=a.terms[i];
  200.         bc=b.terms[i];
  201.         ret.terms[i]=ac-bc;
  202.         }
  203.     return(ret);
  204.     }
  205.  
  206. t_s operator*(t_s &a,t_s &b)
  207.     {
  208.     // Multiply two taylor series
  209.     register int i;
  210.     t_s ret;
  211.     
  212.     set_data(ret,a);
  213.     ret.terms[0]=a.terms[0]*b.terms[0];
  214.     for (i=1;i<a.last_term;i++)
  215.         ret.terms[i]=ats(i,a,b,0);
  216.     return(ret);
  217.     }
  218.  
  219.  
  220. t_s operator/(t_s &a,t_s &b)
  221.     {
  222.     // divide two taylor series
  223.     register int i;
  224.     t_s ret;
  225.     
  226.     set_data(ret,a);
  227.     ret.terms[0]=a.terms[0]/b.terms[0];
  228.     for (i=1;i<a.last_term;i++)
  229.         {
  230.         ret.terms[i]=(a.terms[i]-ats(i,b,ret,1))/b.terms[0];
  231.         }
  232.     return(ret);
  233.     }
  234.  
  235. t_s x(t_s & a)
  236.     {
  237.     // independant variable - x
  238.     register int i;
  239.     t_s ret;
  240.     
  241.     set_data(ret,a);
  242.     ret.terms[0]=a.x0;
  243.     ret.terms[1]=complex(1.0,0.0)*ret.h;
  244.     for (i=2;i<a.last_term;i++)
  245.         ret.terms[i]=complex(0.0,0.0);
  246.     return(ret);
  247.     }
  248.  
  249. t_s constant(t_s & a,complex b)
  250.     {
  251.     // constant - b
  252.     register int i;
  253.     t_s ret;
  254.     
  255.     set_data(ret,a);
  256.     ret.terms[0]=b;
  257.     for (i=1;i<a.last_term;i++)
  258.         ret.terms[i]=complex(0.0,0.0);
  259.     return(ret);
  260.     }
  261.  
  262. t_s expt(t_s & y,t_s & y1)
  263.     {
  264.     // y to power y1
  265.     register int k,ka;
  266.     t_s ret,a1,a2;
  267.     
  268.     set_data(ret,y);
  269.     set_data(a1,y);
  270.     set_data(a2,y);
  271.     ret.terms[0]=expt(y.terms[0],y1.terms[0]);
  272.     a1.terms[0]=lnn(y.terms[0]);
  273.     for (k=1;k<y.last_term;k++)
  274.         {
  275.         ka = k - 1;
  276.         a1.terms[k]=(y.terms[k] - att(ka,y,a1,1))/y.terms[0];
  277.         a2.terms[ka]=ats(k,y,a1,0)*complex(double(k),0.0)/y.h;
  278.         ret.terms[k]=ats(ka,ret,a2,0)*y.h/complex(double(k),0.0);
  279.         }
  280.     return(ret);
  281.     }
  282.  
  283. t_s sqrt(t_s & y)
  284.     {
  285.     // square root
  286.     register int k;
  287.     t_s ret;
  288.     
  289.     set_data(ret,y);
  290.     ret.terms[0]=my_sqrt(y.terms[0]);
  291.     for (k=1;k<y.last_term;k++)
  292.         {
  293.         ret.terms[k]=complex(0.0,0.0);
  294.         ret.terms[k]=(y.terms[k]-ats(k,ret,ret,1))/(ret.terms[0]*
  295.         complex(2.0,0.0));
  296.         }
  297.     return(ret);
  298.     }
  299.     
  300. t_s exp(t_s & y)
  301.     {
  302.     // e to y power
  303.     register int k,ka;
  304.     t_s ret;
  305.     
  306.     set_data(ret,y);
  307.     ret.terms[0]=exp(y.terms[0]);
  308.     for (k=1;k<y.last_term;k++)
  309.         {
  310.         ka = k - 1;
  311.         ret.terms[k]=att(ka,ret,y,0);
  312.         }
  313.     return(ret);
  314.     }
  315. t_s log(t_s & y)
  316.     {
  317.     // ln of y
  318.     register int k,ka;
  319.     t_s ret;
  320.     
  321.     set_data(ret,y);
  322.     ret.terms[0]=lnn(y.terms[0]);
  323.     ret.terms[1]=y.terms[1]/y.terms[0];
  324.     for (k=2;k<y.last_term;k++)
  325.         {
  326.         ka = k - 1;
  327.         ret.terms[k]=(y.terms[k] - att(ka,y,ret,1))/y.terms[0];
  328.         }
  329.     return(ret);
  330.     }
  331.  
  332. t_s sin(t_s & y)
  333.     {
  334.     // sin of y
  335.     register int k,ka;
  336.     t_s f,g;
  337.     
  338.     set_data(f,y);
  339.     set_data(g,y);
  340.     f.terms[0]=sin(y.terms[0]);
  341.     g.terms[0]=cos(y.terms[0]);
  342.     for (k=1;k<y.last_term;k++)
  343.         {
  344.         ka = k - 1;
  345.         f.terms[k]=att(ka,g,y,0);
  346.         g.terms[k]=complex(-1.0,0.0)*att(ka,f,y,0);
  347.         }
  348.     return(f);
  349.     }
  350.     
  351. t_s cos(t_s & y)
  352.     {
  353.     // cos of y
  354.     register int k,ka;
  355.     t_s f,g;
  356.     
  357.     set_data(f,y);
  358.     set_data(g,y);
  359.     f.terms[0]=sin(y.terms[0]);
  360.     g.terms[0]=cos(y.terms[0]);
  361.     for (k=1;k<y.last_term;k++)
  362.         {
  363.         ka = k - 1;
  364.         f.terms[k]=att(ka,g,y,0);
  365.         g.terms[k]=complex(-1.0,0.0)*att(ka,f,y,0);
  366.         }
  367.     return(g);
  368.     }
  369.     
  370. t_s tan(t_s & y)
  371.     {
  372.     // cos of y
  373.     register int k,ka;
  374.     t_s f,a1,a2;
  375.     
  376.     set_data(f,y);
  377.     set_data(a1,y);
  378.     set_data(a2,y);
  379.     a1.terms[0]=sin(y.terms[0]);
  380.     a2.terms[0]=cos(y.terms[0]);
  381.     f.terms[0]=a1.terms[0]/a2.terms[0];
  382.     for (k=1;k<y.last_term;k++)
  383.         {
  384.         ka = k - 1;
  385.         a1.terms[k]=att(ka,a2,y,0);
  386.         a2.terms[k]=complex(-1.0,0.0)*att(ka,a1,y,0);
  387.         f.terms[k]=(a1.terms[k] - ats(k,a2,f,1))/a2.terms[0];
  388.         }
  389.     return(f);
  390.     }
  391.     
  392. t_s sinh(t_s & y)
  393.     {
  394.     // sinh of y
  395.     register int k,ka;
  396.     t_s f,g;
  397.     
  398.     set_data(f,y);
  399.     set_data(g,y);
  400.     f.terms[0]=sinh(y.terms[0]);
  401.     g.terms[0]=cosh(y.terms[0]);
  402.     for (k=1;k<y.last_term;k++)
  403.         {
  404.         ka = k - 1;
  405.         f.terms[k]=att(ka,g,y,0);
  406.         g.terms[k]=att(ka,f,y,0);
  407.         }
  408.     return(f);
  409.     }
  410.     
  411. t_s cosh(t_s & y)
  412.     {
  413.     // cosh of y
  414.     register int k,ka;
  415.     t_s f,g;
  416.     
  417.     set_data(f,y);
  418.     set_data(g,y);
  419.     f.terms[0]=sinh(y.terms[0]);
  420.     g.terms[0]=cosh(y.terms[0]);
  421.     for (k=1;k<y.last_term;k++)
  422.         {
  423.         ka = k - 1;
  424.         f.terms[k]=att(ka,g,y,0);
  425.         g.terms[k]=att(ka,f,y,0);
  426.         }
  427.     return(g);
  428.     }
  429.     
  430. t_s tanh(t_s & y)
  431.     {
  432.     // tanh of y
  433.     register int k,ka;
  434.     t_s f,a1,a2;
  435.     
  436.     set_data(f,y);
  437.     set_data(a1,y);
  438.     set_data(a2,y);
  439.     a1.terms[0]=sinh(y.terms[0]);
  440.     a2.terms[0]=cosh(y.terms[0]);
  441.     f.terms[0]=a1.terms[0]/a2.terms[0];
  442.     for (k=1;k<y.last_term;k++)
  443.         {
  444.         ka = k - 1;
  445.         a1.terms[k]=att(ka,a2,y,0);
  446.         a2.terms[k]=att(ka,a1,y,0);
  447.         f.terms[k]=(a1.terms[k] - ats(k,a2,f,1))/a2.terms[0];
  448.         }
  449.     return(f);
  450.     }
  451.     
  452. t_s asin(t_s & y)
  453.     {
  454.     // asin of y
  455.     register int k,ka;
  456.     t_s f,a1;
  457.     
  458.     set_data(f,y);
  459.     set_data(a1,y);
  460.     f.terms[0]=asin(y.terms[0]);
  461.     a1.terms[0]=cos(f.terms[0]);
  462.     for (k=1;k<y.last_term;k++)
  463.         {
  464.         ka = k - 1;
  465.         f.terms[k]=(y.terms[k] - att(ka,a1,f,1))/a1.terms[0];
  466.         a1.terms[k]=complex(-1.0,0.0)*att(ka,y,f,0);
  467.         }
  468.     return(f);
  469.     }
  470.     
  471. t_s acos(t_s & y)
  472.     {
  473.     // acos of y
  474.     register int k,ka;
  475.     t_s f,a1;
  476.     
  477.     set_data(f,y);
  478.     set_data(a1,y);
  479.     f.terms[0]=acos(y.terms[0]);
  480.     a1.terms[0]=sin(f.terms[0]);
  481.     for (k=1;k<y.last_term;k++)
  482.         {
  483.         ka = k - 1;
  484.         f.terms[k]=complex(-1.0,0.0)*
  485.         (y.terms[k] + att(ka,a1,f,1))/a1.terms[0];
  486.         a1.terms[k]=att(ka,y,f,0);
  487.         }
  488.     return(f);
  489.     }
  490.     
  491. t_s atan(t_s & y)
  492.     {
  493.     // atan of y
  494.     register int k,ka;
  495.     t_s f,a1,a2;
  496.     
  497.     set_data(f,y);
  498.     set_data(a1,y);
  499.     set_data(a2,y);
  500.     f.terms[0]=atan(y.terms[0]);
  501.     a1.terms[0]=sin(f.terms[0]);
  502.     a2.terms[0]=cos(f.terms[0]);
  503.     for (k=1;k<y.last_term;k++)
  504.         {
  505.         ka = k - 1;
  506.         f.terms[k]=(ats(k,y,a2,1)-y.terms[0]*att(ka,a1,f,1) -
  507.             att(ka,a2,f,1))/(a2.terms[0]+y.terms[0]*a1.terms[0]);
  508.         a1.terms[k]=att(ka,a2,f,0);
  509.         a2.terms[k]=complex(-1.0,0.0)*att(ka,a1,f,0);
  510.         }
  511.     return(f);
  512.     }
  513.  
  514. t_s deriv(t_s &a) 
  515.     {
  516.     // derivative of a taylor series
  517.     register int i;
  518.     t_s ret;
  519.         
  520.     set_data(ret,a);
  521.     
  522.     for (i=a.first_term;i<a.last_term-1;i++)
  523.         {
  524.         if (i > 0)
  525.             ret.terms[i]=a.terms[i+1]*complex(double(i),0.0)/a.h;
  526.         else 
  527.             ret.terms[i]=a.terms[i+1]/a.h;
  528.         }
  529.     ret.terms[a.last_term]=complex(0.0,0.0); // unknown but small
  530.     return(ret);
  531.     }
  532. void diff_eq(t_s &ret2,t_s &y)
  533.     {
  534.     // solve 1st order ordinary differential equation provided in "equation".
  535.     register int k;
  536.     t_s work;
  537.     t_s work2;
  538.     
  539.     work=y;
  540.     work2=y;
  541.     for (k=work.last_term;k<MAX_TERMS;k++)
  542.         {
  543.         // get derivative t_s (work) from function t_s work2
  544.         equation(work,work2); 
  545.         set_term_from_deriv(k,work,work2);
  546.         }
  547.     ret2 = work2;
  548.     }
  549. void higher_diff_eq(t_s &ret2,t_s &y,int n)
  550.     {
  551.     // solve 1st order ordinary differential equation provided in "equation".
  552.     register int k;
  553.     t_s work;
  554.     t_s work2;
  555.     
  556.     work=y;
  557.     work2=y;
  558.     for (k=work.last_term;k<MAX_TERMS;k++)
  559.         {
  560.         // get derivative t_s (work) from function t_s work2
  561.         equation(work,work2); 
  562.         higher_set_term_from_deriv(k,work,work2,n);
  563.         }
  564.     ret2 = work2;
  565.     }
  566. void set_term_from_deriv(int i,t_s &dy,t_s &y)
  567.     {
  568.     if (i > 0) 
  569.         y.terms[i]=dy.terms[i-1]*y.h/complex(double(i),0.0);
  570.     else 
  571.         y.terms[i]=dy.terms[i-1]*y.h;
  572.     y.last_term++;
  573.     }
  574. void higher_set_term_from_deriv(int i,t_s &dy,t_s &y,int n)
  575.     {
  576.     complex adj1,adj2;
  577.     
  578.     adj1 = complex(1.0,0.0);
  579.     adj2 = complex(1.0,0.0);
  580.     
  581.     for (int a=0;a<n;a++)
  582.         adj1 *= dy.h;
  583.     for (int b=0;b<n;b++)
  584.         {
  585.         if (i+b > 0)
  586.             {
  587.             adj2 *= complex(double(i+b),0.0);
  588.             }
  589.         }
  590.     y.terms[i]=dy.terms[i-n]*adj1/adj2;
  591.     y.last_term++;
  592.     }
  593. void jump(t_s &y)
  594.     {
  595.     t_s y2;
  596.     
  597.     y2     = y;
  598.     y2.terms[0] = complex(0.0,0.0);
  599.     for (int i = 0;i < y.last_term;i++)
  600.         {
  601.         y2.terms[0] += y.terms[i];
  602.         }
  603.     y2.x0 = y.x0 + y.h;
  604.     y2.last_term = 1;
  605.     y = y2;
  606.     }
  607.  
  608. void leap(t_s &y,int n)
  609.     {
  610.     t_s y2;
  611.     complex adj1,adj2;
  612.         
  613.     y2     = y;
  614.  
  615.     adj2 = complex(1.0,0.0); 
  616.     for (int j = 0;j < n; j++)
  617.         {
  618.         y2.terms[j] = complex(0.0,0.0);
  619.         for (int i = 0;i < y.last_term-j;i++)
  620.             {
  621.             adj1 = complex(1.0,0.0);
  622.             for (int k = 1 ; k <= j; k++)
  623.                 {
  624.                 adj1 *= complex(double(i+k),0.0);
  625.                 }
  626.             y2.terms[j] += y.terms[i+j]*adj1;
  627.             }
  628.         if (j > 0)
  629.             {
  630.             adj2 *= j;
  631.             y2.terms[j] /= adj2;
  632.             }
  633.         }
  634.     y2.x0 = y.x0 + y.h;
  635.     y2.last_term = n;
  636.     y = y2;
  637.     }
  638.  
  639. int get_last_term(t_s &y)
  640.     {
  641.     return(y.last_term);
  642.     }
  643.     
  644. double get_x0(t_s &y)
  645.     {
  646.     return(real(y.x0));
  647.     }
  648. double get_y0(t_s &y)
  649.     {
  650.     return(real(y.terms[0]));
  651.     }
  652.  
  653. void adjust_h(t_s &y,complex new_h)
  654.     {
  655.     t_s y2;
  656.     complex t,adj;
  657.     
  658.     t = new_h / y.h;
  659.     adj = t;
  660.     y2     = y;
  661.     for (int i = 1;i < y.last_term;i++)
  662.         {
  663.         y2.terms[i] = y.terms[i]*adj;
  664.         adj *= t;
  665.         }
  666.     y2.h = new_h;
  667.     y = y2;
  668.     }
  669.