dsolve.cpp

   1 #include "stdafx.h"
   2
   3 #include "defs.h"
   4
   5 //      q(x)y' + p(x)*y = g(x)
   6 //
   7 //      u(x) = exp(integral(p))
   8 //
   9 //      y = (integral(u*g) + c) / u(x)
  10
  11
  12 #define f p1
  13 #define y p2
  14 #define x p3
  15
  16 #define p p4
  17 #define g p5
  18 #define q p6
  19
  20 #define mu p7
  21
  22 void
  23 dsolve(void)
  24 {
  25         int n;
  26
  27         save();
  28
  29         x = pop();
  30         y = pop();
  31         f = pop();
  32
  33         push(f);
  34         push(y);
  35         push(x);
  36
  37         n = distilly();
  38
  39         if (n != 3)
  40                 stop("error in dsolve");
  41
  42         q=pop();
  43
  44         p = pop();
  45
  46         negate();
  47         g = pop();
  48
  49 /*      print(g);
  50         print(p);
  51         print(p);
  52 */
  53         push(p);
  54         push(q);
  55         divide();
  56         push(x);
  57         integral();
  58         exponential();
  59         mu = pop();
  60
  61         push(mu);
  62         push(g);
  63         push(q);
  64         divide();
  65         multiply();
  66         push(x);
  67         integral();
  68         scan("C");
  69         add();
  70         push(mu);
  71         divide();
  72
  73         restore();
  74 }
  75
  76 //      n       p1              p2      p3      p4      p5      stack
  77
  78 //      1       4y'+3xy+2x+1    y       x       1       2x+1    2x+1
  79
  80 //      2       4y'+3xy         y'      x       y       3xy     3x
  81
  82 //      3       4y'             y''     x       y'      4y'     4
  83
  84 int distilly()
  85 {
  86         int n = 0;
  87         save();
  88         p4 = one;
  89         p3 = pop();
  90         p2 = pop();
  91         p1 = pop();
  92         while (1) {
  93                 n++;
  94                 push(p1);
  95                 push(p2);
  96                 push(zero);
  97                 subst();
  98                 eval();
  99                 p5 = pop();
 100                 push(p5);
 101                 push(p4);
 102                 divide();
 103                 push(p1);
 104                 push(p5);
 105                 subtract();
 106                 p1 = pop();
 107                 if (equal(p1, zero))
 108                         break;
 109                 p4 = p2;
 110                 push(p2);
 111                 push(p3);
 112                 derivative();
 113                 p2 = pop();
 114         }
 115         restore();
 116         return n;
 117 }