share/mnewton/rtest_mnewton.mac

   1 ( kill (a, b, c, d, u, v, x, x1, x2, y, z),
   2   approx_equal(x,y,epsilon) := block(
   3     [vx:setify(listofvars(x)),vy:setify(listofvars(y)),prederror:true],
   4     if vx#vy then return(false),
   5     x : sort(flatten(listify(x))), y : sort(flatten(listify(y))),
   6     every (lambda ([u, v], is (abs (rectform (rhs(u) - rhs(v))) < epsilon)), x, y)),
   7   if errcatch (fundef (mnewton)) = [] then load("mnewton"),
   8   0);
   9 0$
  10
  11 approx_equal(
  12   mnewton(x-1.0,x,0),
  13   [[x = 1.0]],
  14   newtonepsilon
  15   )$
  16 true $
  17
  18 approx_equal(
  19   mnewton(sin(x),x,1.0),
  20   [[x = 0.0]],
  21   newtonepsilon
  22   )$
  23 true $
  24
  25
  26 /* from manual */
  27 approx_equal(
  28   mnewton([x1+3*log(x1)-x2^2, 2*x1^2-x1*x2-5*x1+1], [x1, x2], [5, 5]),
  29   [[x1 = 3.756834008012769, x2 = 2.779849592817897]],
  30   newtonepsilon
  31   )$
  32 true $
  33
  34 /* from manual */
  35 approx_equal(
  36   mnewton([2*a^a-5],[a],[1]),
  37   [[a = 1.70927556786144]],
  38   newtonepsilon
  39   )$
  40 true $
  41
  42 /* from manual */
  43 approx_equal(
  44   mnewton([2*3^u-v/u-5, u+2^v-4], [u, v], [2, 2]),
  45   [[u = 1.066618389595407, v = 1.552564766841786]],
  46   newtonepsilon
  47   )$
  48 true $
  49
  50 /* same as preceding example, with option argument DF */
  51 approx_equal (
  52   mnewton(eqs: [2*3^u-v/u-5, u+2^v-4], [u, v], [2, 2], jacobian (eqs, [u, v])),
  53   [[u = 1.066618389595407, v = 1.552564766841786]],
  54   newtonepsilon) $
  55 true $
  56
  57 /* see https://sourceforge.net/p/maxima/bugs/1199/ */
  58 approx_equal(
  59   mnewton([2*(x1-1.1234),3*(x2-2.2345)], [x1, x2], [1,2]),
  60   [[x1=1.1234,x2=2.2345]],
  61   newtonepsilon
  62   )$
  63 true $
  64
  65
  66 block([fpprec: 16, newtonepsilon:bfloat(newtonepsilon)],
  67   approx_equal(
  68     mnewton(cos(x)-%i,x,1b0),
  69     [[x = 1.570796326794897b0 - 8.81373587019543b-1*%i]],
  70     newtonepsilon
  71     ))$
  72 true $
  73
  74 /* see http://www.math.utexas.edu/pipermail/maxima/2010/021238.html and
  75 https://sourceforge.net/p/maxima/bugs/1967/ */
  76 block([f:diff(acos(x^2+8.363+267)/(sqrt(x^2-4.29*x+1042)*(x^2+12.8102*x+1177)),x)],
  77   approx_equal(
  78     mnewton(f,x,-5.0),
  79     [[x = 5.69311809542034e-9*%i - 5.626466372932345]],
  80     newtonepsilon
  81     ))$
  82 true $
  83
  84 block([fpprec:100, newtonepsilon:bfloat(newtonepsilon)*1b-8,f:diff(acos(x^2+8.363+267)/(sqrt(x^2-4.29*x+1042)*(x^2+12.8102*x+1177)),x)],
  85   approx_equal(
  86     mnewton(f,x,-5.0),
  87     [[x = - 5.626466372928149b0]],
  88     newtonepsilon
  89     ))$
  90 true $
  91
  92 block([fpprec:100,newtonepsilon:bfloat(10^(-fpprec/2)),f:diff(acos(x^2+8.363+267)/(sqrt(x^2-4.29*x+1042)*(x^2+12.8102*x+1177)),x)],
  93   approx_equal(
  94     mnewton(f,x,-5.0),
  95     [[x = - 5.62646637292814898636208310765236683228303380080106742342217423415231803837099458924175436056796413b0]],
  96     newtonepsilon
  97     ))$
  98 true $
  99
 100 /* see https://sourceforge.net/p/maxima/bugs/1133/ */
 101 approx_equal(
 102   mnewton([x^y - %i,x+y-2],[x,y],[1,1.5]),
 103   [[x=0.563519172421*%i+0.451201761049,y=1.5487982389505-0.563519172421*%i]],
 104   newtonepsilon
 105   )$
 106 true $
 107
 108 /* see https://sourceforge.net/p/maxima/bugs/1118/ */
 109 approx_equal(
 110   mnewton([x[1] + x[2]-1, x[1] - x[2] - 10],[x[1],x[2]],[1,2]),
 111   [[x[1] = 5.5,x[2] = -4.5]],
 112   newtonepsilon
 113   )$
 114 true $
 115
 116 /* A test for a command that failed in 2010, see
 117  https://def.fe.up.pt/pipermail/maxima-discuss/2010/033075.html and
 118  https://def.fe.up.pt/pipermail/maxima-discuss/2010/033080.html
 119
 120  As the maxima-discuss thread "COERCE-FLOAT-FUN in complexfield"
 121  in Mid-Jan 2019 shows the remedy for this made the following
 122  command slow and space-hungry in sbcl:
 123
 124      w[i,j] := random (1.0) + %i * random (1.0);
 125      showtime : true$
 126      M : genmatrix (w, 100, 100)$
 127      lu_factor (M, complexfield)$
 128      lu_factor (M, generalring)$
 129
 130  The reason for that was using coerce-float-fun for each matrix
 131  entry which created lisp functions sbcl automatically compiled.
 132 */
 133 approx_equal(
 134   mnewton(diff(acos(x^2+8.363+267)/(sqrt(x^2-4.29*x+1042)*(x^2+12.8102*x+1177)),x),x,5),
 135   [[x=-1.012828847267113*10^-12*%i-5.626466372928149]],
 136   newtonepsilon
 137 );
 138 true $
 139
 140 /* maxima-discuss 2020-06-21: "mnewton error message" */
 141
 142 block ([s: make_string_output_stream ()],
 143        with_stdout (s, mnewton([x^2+y^2+z^2,x^2-y^2+z^2],[x,y],[1,1])),
 144        [%%, get_output_stream_string (s)]);
 145 [false, "mnewton: extra variables in equations; found: z
 146 "];
 147
 148 /* verify input representation options */
 149
 150 approx_equal (mnewton ([x^2 - y^2 - 5, x*y - 2], [x, y], [1, 1]), [[x = 2.387794404616198, y = 0.8375930507808815]], newtonepsilon);
 151 true;
 152
 153 approx_equal (mnewton ([x^2 - y^2 = 5, x*y = 2], [x, y], [1, 1]), [[x = 2.387794404616198, y = 0.8375930507808815]], newtonepsilon);
 154 true;
 155
 156 approx_equal (mnewton (- [x^2 - y^2 - 5, x*y - 2], [x, y], [1, 1]), [[x = 2.387794404616198, y = 0.8375930507808815]], newtonepsilon);
 157 true;
 158
 159 approx_equal (mnewton ([5 = x^2 - y^2, 2 = x*y], [x, y], [1, 1]), [[x = 2.387794404616198, y = 0.8375930507808815]], newtonepsilon);
 160 true;
 161
 162 approx_equal (mnewton ([x[1]^2 - x[2]^2 - 5, x[1]*x[2] - 2], [x[1], x[2]], [1, 1]), [[x[1] = 2.387794404616198, x[2] = 0.8375930507808815]], newtonepsilon);
 163 true;
 164
 165 approx_equal (mnewton ([z^3 - 1], [z], [1 + %i]), [[z = 1.0]], newtonepsilon);
 166 true;
 167
 168 approx_equal (mnewton ([z^3 = 1], [z], [1 + %i]), [[z = 1.0]], newtonepsilon);
 169 true;
 170
 171 approx_equal (mnewton (z^3 - 1, z, 1 + %i), [[z = 1.0]], newtonepsilon);
 172 true;
 173
 174 approx_equal (mnewton (z^3 = 1, z, 1 + %i), [[z = 1.0]], newtonepsilon);
 175 true;
 176
 177
 178 /* additional test cases */
 179
 180 block ([eqs: [sin(%pi*(b + a)) - a^2, b^2 - %e*a], eqs1, aa: %pi, bb: 5],
 181        eqs1: map ("=", eqs, subst ([a = aa, b = bb], eqs)),
 182        approx_equal (mnewton (eqs1, [a, b], [1, 1]), [[a = aa, b = bb]], newtonepsilon));
 183 true;
 184
 185 block ([eqs: [sqrt(b^3 + a) - a, exp(b) - a/%pi], eqs1, aa: 0, bb: 1],
 186        eqs1: map ("=", eqs, subst ([a = aa, b = bb], eqs)),
 187        approx_equal (mnewton (eqs1, [a, b], [1, 1]), [[a = aa, b = bb]], newtonepsilon));
 188 true;
 189
 190 block ([eqs: [sqrt(b^3 + a) - a*c, exp(b) - a/%pi - c, c^2 + a - b], eqs1, aa: 0, bb: 1, cc: 2],
 191        eqs1: map ("=", eqs, subst ([a = aa, b = bb, c = cc], eqs)),
 192        approx_equal (mnewton (eqs1, [a, b, c], [1, 1, 1]), [[a = aa, b = bb, c = cc]], newtonepsilon));
 193 true;
 194
 195 block ([eqs: [sin(%pi*(b + a)) + a^2*%i, b^2 + %e*a], eqs1, aa: %pi, bb: - 5],
 196        eqs1: map ("=", eqs, subst ([a = aa, b = bb], eqs)),
 197        approx_equal (mnewton (eqs1, [a, b], [%i, %i]), [[a = aa, b = bb]], newtonepsilon));
 198 true;
 199
 200 block ([eqs: [sqrt(b^3 + a) - a*c + d, exp(b) - a/%pi - c*d, c^2 + a - b, b*log(d) - a*c], eqs1,
 201         aa: 0.9316994624088384 - 0.4532672590238549*%i,
 202         bb: 0.03715752272566857*%i + 0.07362074010484573,
 203         cc: 1.778327483615696 - 0.1432737285301378*%i,
 204         dd: 2.99576612219564*%i - 0.1070840677802114],
 205        eqs1: map ("=", eqs, subst ([a = aa, b = bb, c = cc, d = dd], eqs)),
 206        approx_equal (mnewton (eqs1, [a, b, c, d], [1, %i, 1, %i]), [[a = aa, b = bb, c = cc, d = dd]], newtonepsilon));
 207 true;
 208