share/contrib/diffequations/tests/rtest_ode1_abel.mac

   1 (load("contrib_ode"),0);
   2 0$
   3
   4 /* Test cases for ode1_abel.mac functions */
   5
   6
   7 /* Print ode number*/
   8 (pn_(n_):=print("ODE ",n_),kn_(n_):=print("Kamke ODE ",n_),true);
   9 true;
  10
  11 /* Murphy 78 Abel */
  12 (pn_(78),ans:ode1_abel(eqn:'diff(y,x)+(a*x+y)*y^2=0,y,x));
  13 [];
  14
  15 /* Murphy 79 Abel */
  16 (pn_(79),ans:ode1_abel(eqn:'diff(y,x)+(a*exp(x)+y)*y^2=0,y,x));
  17 [];
  18
  19 /* about the simplest Abel equation */
  20 ans:ode1_abel(eqn:'diff(y,x)=(x+c)*y^3+y^2,y,x);
  21 [-((sqrt(3)*log((x+c)^2*y^2+(x+c)*y+1)
  22     +2*atan((2*(x+c)*y+1)/sqrt(3))-2*sqrt(3)*log((x+c)*y))
  23  /(2*sqrt(3))) = log(x+c)+%c];
  24 [method,ode_check(eqn,ans[1])];
  25 [abel,0];
  26
  27 ans:ode1_abel(eqn:'diff(y,x)=2*x^3*y^3+2*x*y^2,y,x);
  28 [-((sqrt(3)*log(x^4*y^2+x^2*y+1)
  29     +2*atan((2*x^2*y+1)/sqrt(3))
  30     -2*sqrt(3)*log(x^2*y))
  31  /(4*sqrt(3))) = log(x)+%c];
  32 [method,ode_check(eqn,ans[1])];
  33 [abel,0];
  34
  35 /* and now some cases with relative invariant = 0.  All generated by solving
  36    for f[0] for "random" f[1],f[2],f[3]
  37 */
  38 ans:ode1_abel(eqn:'diff(y,x,1)=x*y^3+x*y^2+y+(9*x^2-2*x^3)/(27*x^2),y,x);
  39 [y = %e^(x-x^2/6)/sqrt(%c-sqrt(3)*%e^3*%i
  40  *((sqrt(3)*%i*gamma_incomplete(1,(3*(2-(2*x)/3)^2)/4)*(2-(2*x)/3)^2)/((2*x)/3-2)^2
  41    -(3*%i*gamma_incomplete(1/2,(3*(2-(2*x)/3)^2)/4)*(2-(2*x)/3))/abs((2*x)/3-2)))
  42    -1/3]$
  43 ode_check(eqn,ans[1]);
  44 0;
  45
  46 ans:ode1_abel(eqn:'diff(y,x,1)=y^3/x+y,y,x);
  47 [y = %e^x/sqrt(%c+2*gamma_incomplete(0,-2*x))];
  48 ode_check(eqn,ans[1]);
  49 0;
  50
  51 ans:ode1_abel(eqn:'diff(y,x,1)=y^3+x^2*y,y,x);
  52 [y=%e^(x^3/3)/sqrt(%c-2*3^(1/3)*gamma_incomplete(1/3,-(2*x^3/3))/(3*2^(1/3)))];
  53 ode_check(eqn,ans[1]);
  54 0;
  55
  56 ans:ode1_abel(eqn:'diff(y,x,1)=x^2*y^3+x*y,y,x);
  57 [y=%e^(x^2/2)/sqrt(%c-2*(sqrt(%pi)*%i*erf(%i*x)/4+x*%e^x^2/2))];
  58 ode_check(eqn,ans[1]);
  59 0;
  60
  61 /* Cheb-Terrab and Roche (31) */
  62 ans:ode1_abel(eqn:'diff(y,x)=y^3/x+y^2,y,x);
  63 [];
  64
  65 /* Cheb-Terrab and Roche (33) */
  66 ans:ode1_abel(eqn:'diff(y,x)=(x+1/x+1/x^3)*y^3/x+y^2,y,x);
  67 [];
  68
  69 /* According to Cheb-Terrab and Roche, the following
  70    equations from kamke have constant invariant */
  71
  72 /* Kamke 38 - Chini invariant constant - 1/(8*a*b^2)
  73
  74    One integrable case is a:-1/6  b:1/3
  75 */
  76
  77 (kn_(38),a:-1/6,b:1/3,ans:ode1_abel(eqn:'diff(y,x) - a*y^3 - b*x^(-3/2),y,x));
  78 [-log(3*sqrt(x)*y+3)/3+log(3*sqrt(x)*y/2-3)/3+3/(3*sqrt(x)*y+3)+3*log(x)/2+log(1/(54*x^3))/3=%c];
  79 [method,ode_check(eqn,ans[1])];
  80 [abel,0];
  81 kill(a,b);
  82 done;
  83
  84 /* Kamke 41 f0=f1=0 */
  85 assume(b^2>0,(b^2+4*a)>0);
  86 [b^2>0,(b^2+4*a)>0];
  87 (kn_(41),ans:ode1_abel(eqn:'diff(y,x) + a*x*y^3 + b*y^2 ,y,x));
  88 [-((b*log(-((-2*a*x*y+sqrt(b^2+4*a)-b)/(2*a*x*y+sqrt(b^2+4*a)+b)))
  89    +sqrt(b^2+4*a)*(log(a^2*x^2*y^2+a*b*x*y-a)-2*log(a*x*y/b)))
  90  /(2*sqrt(b^2+4*a))) = log(x)+%c];
  91 [method,ode_check(eqn,ans[1])];
  92 [abel,0];
  93 forget(b^2>0,(b^2+4*a)>0);
  94 [b^2>0,(b^2+4*a)>0];
  95
  96 /* Kamke 46  relative invariant = 0 */
  97 (kn_(46),a:1,ans:ode1_abel(eqn:'diff(y,x)- x^a*y^3+3*y^2-x^(-a)*y-x^(-2*a)+ a*x^(-a-1),y,x));
  98 [y=1/x+1/(sqrt(1/x^2+%c)*x^2)];
  99 ode_check(eqn,ans[1]);
 100 0;
 101 kill(a);
 102 done;
 103 assume(a>1);
 104 [a>1];
 105 (kn_(46),ans:ode1_abel(eqn:'diff(y,x)- x^a*y^3+3*y^2-x^(-a)*y-x^(-2*a)+ a*x^(-a-1),y,x));
 106 [y = %e^-(2*x^(1-a)/(1-a))/sqrt(2*(a-1)^((a+1)/(1-a))
 107           *gamma_incomplete((a+1)/(1-a),-4*x^(1-a)/(a-1))*x^(a+1)
 108   /((1-a)*4^((a+1)/(1-a))*(-x^(1-a))^((a+1)/(1-a)))+%c)+1/x^a];
 109 ode_check(eqn,ans[1]);
 110 0;
 111 forget(a>1);
 112 [a>1];
 113
 114 /* Kamke 49 */
 115 /*(kn_(49),ans:ode1_abel(eqn:'diff(y,x)+a*diff(phi(x),x)*y^3+6*a*phi(x)*y^2+(2*a+1)*y*diff(phi(x),x,2)/diff(phi(x),x) +2*(a+1) ,y,x));
 116 false;*/
 117
 118 /* Kamke 51 */
 119 /*(kn_(51),ans:ode1_abel(eqn:'diff(y,x)-(y-f(x))*(y-g(x))*(y-(a*f(x)+b*g(x))/(a+b))*h(x)-'diff(f(x),x)*(y-g(x))/(f(x)-g(x))-'diff(g(x),x)*(y-f(x))/(g(x)-f(x)),y,x));
 120 false;*/
 121
 122 /* Kamke 188 */
 123 /* Some choices that are integrable include
 124     (n:3, b:1, a:n+b);   => K = -27/4
 125     (n:7, b:2, a:n+b);   => K = -343/36
 126 */
 127 (kn_(188),n:3,b:1,a:n+b,ans:ode1_abel(eqn:x^(2*n+1)*'diff(y,x)-a*y^3-b*x^(3*n),y,x));
 128 [log(-3*y/x^3-3)/3-log(3-6*y/x^3)/3+3/(3-6*y/x^3)-2*log(x)+log(-4/x^3)/3=%c];
 129 [method,ode_check(eqn,ans[1])];
 130 [abel,0];
 131 kill(a,b,n);
 132 done;
 133
 134 /* Non-constant invariant */
 135 (kn_(36),ans:ode1_abel(eqn:'diff(y,x)+y^3+a*x*y^2,y,x));
 136 [];
 137
 138 (kn_(37),ans:ode1_abel(eqn:'diff(y,x)-y^3-a*exp(x)*y^2,y,x));
 139 [];
 140
 141 (kn_(40),ans:ode1_abel(eqn:'diff(y,x)+3*a*y^3+6*a*x*y^2,y,x));
 142 [];
 143
 144 (kn_(42),ans:ode1_abel(eqn:'diff(y,x)-x*(x+2)*y^3-(x+3)*y^2,y,x));
 145 [];
 146
 147 (kn_(43),ans:ode1_abel(eqn:'diff(y,x)+(3*a*x^2+4*a^2*x+b)*y^3+3*x*y^2,y,x));
 148 [];
 149
 150 (kn_(45),ans:ode1_abel(eqn:'diff(y,x)+2*(a^2*x^3-b^2*x)*y^3+3*b*y^2,y,x));
 151 [];
 152
 153 /* Abel equations of second kind */
 154 (kn_(213),ans:ode1_abel(eqn:(y+1)*'diff(y,x)=y+x,y,x));
 155 [-((log((2*(x-1)/(y+1)-sqrt(5)+1)/(2*(x-1)/(y+1)+sqrt(5)+1))+sqrt(5)*log((x-1)/(y+1)+(x-1)^2/(y+1)^2-1)-2*sqrt(5)*log((x-1)/(y+1)))/(2*sqrt(5)))=log(x-1)+%c];
 156 [method,ode_check(eqn,ans[1])];
 157 [abel2,0];
 158
 159 (kn_(214),ans:ode1_abel(eqn:(y+x-1)*'diff(y,x)-y+2*x+3=0,y,x));
 160 [-((2*atan((1-(3*x+2)/(y+x-1))/sqrt(2))
 161     +sqrt(2)*log(-(2*(3*x+2)/(y+x-1))+(3*x+2)^2/(y+x-1)^2+3)
 162     -2*sqrt(2)*log(-((3*x+2)/(2*(y+x-1)))))
 163  /(6*sqrt(2))) = log(3*x+2)/3+%c];
 164 [method,ode_check(eqn,ans[1])];
 165 [abel2,0];
 166
 167 (kn_(215),ans:ode1_abel(eqn:(y+2*x-2)*'diff(y,x)-y+x+1=0,y,x));
 168 [-(log(-((3*x-1)/(y+2*x-2))+(3*x-1)^2/(3*(y+2*x-2)^2)+1)
 169         +2*sqrt(3)*atan((3-(2*(3*x-1))/(y+2*x-2))/sqrt(3))
 170         -(2*log(-((3*x-1)/(3*(y+2*x-2))))))
 171         /6
 172          = log(3*x-1)/3+%c]$
 173 [method,ode_check(eqn,ans[1])];
 174 [abel2,0];
 175
 176 (kn_(216),ans:ode1_abel(eqn:(y-2*x+1)*'diff(y,x)+y+x=0,y,x));
 177 [-(((2*sqrt(3)*atan(((2*(3*x-1))/(y-2*x+1)+3)/sqrt(3))
 178         +log((3*x-1)/(y-2*x+1)+(3*x-1)^2/(3*(y-2*x+1)^2)+1)
 179         -(2*log((3*x-1)/(3*(y-2*x+1)))))
 180         /6))
 181          = log(3*x-1)/3+%c]$
 182 [method,ode_check(eqn,ans[1])];
 183 [abel2,0];
 184
 185 (kn_(217),ans:ode1_abel(eqn:(y-x^2)*'diff(y,x)=x,y,x));
 186 [(y-x^2)*(log(-1/(2*(y-x^2)))/(2*(y-x^2))-log(1-1/(2*(y-x^2)))/(2*(y-x^2))-1)/2= x^2/2+%c];
 187 [method,ode_check(eqn,ans[1])];
 188 [abel2,0];
 189
 190 (kn_(218),ans:ode1_abel(eqn:(y-x^2)*'diff(y,x)+4*x*y=0,y,x));
 191 [(log(2*x^2/(y-x^2)+2)-2*log(2*x^2/(y-x^2)+1)+log(2*x^2/(3*(y-x^2))))/2=log(x)+%c];
 192 [method,ode_check(eqn,ans[1])];
 193 [abel2,0];
 194
 195 /* Why does this fail? */
 196 /* (kn_(220),ans:ode1_abel(eqn:2*y*'diff(y,x)-x*y^2-x^3=0,y,x));
 197 [%e^-(x^2/2)*(y^2+x^2+2) = %c];
 198 [method,ode_check(eqn,ans[1])];
 199 [abel2,0]; */
 200
 201 (kn_(221),ans:ode1_abel(eqn:(2*y+x+1)*'diff(y,x)-(2*y+x-1)=0,y,x));
 202 [-((2*y+x+1)
 203     * (-(4*(2*(x+1)+2*(1-x))*log(1-(2*(x+1)+2*(1-x))/(3*(2*y+x+1)))/(3*(2*y+x+1)))
 204         +4*(2*(x+1)+2*(1-x))*log(-((2*(x+1)+2*(1-x))/(3*(2*y+x+1))))/(3*(2*y+x+1))-4 )
 205  / (3*(2*(x+1)+2*(1-x)))) = x+%c];
 206 [method,ode_check(eqn,ans[1])];
 207 [abel2,0];