archive/books/schelter/xplot.bk

   1 \x06\x01\x19\x16\x05
   2 ((face (xplot-eval 129 577 580 1139 1163 1612 1650 1830 1837 1968 1978 2166 2271 2496 2544 3036 3198 4078 4091 4505 4517 4647 4660 4740 4753 4993 5006 5441 5454 5757 5767 5899 5909 6041 6050 6709 6724 7209 7225 7695 7705 7901 7915 8402 8417 8864 8878 9222 9231 9568 9578 10379)) (book-command-arg))\f
   3       This file contains a number of 2 and 3 dimensional Plots
   4
   5 Some are from differential equations, others are just shapes.
   6
   7
   8 # Solve the differential equation system
   9 #
  10 #  x' = y
  11 #  y' = -sin(x) - y
  12 # and then plot  [x[t],y[t]] as a function of t.
  13 # in this case there is no z coordinate.
  14
  15 function definitions {
  16           xdot1(x,y)=           y;
  17           xdot2(x,y)=-sin(x) - .1*y;};
  18    sol = object{
  19         ode { [xdot1(x,y),xdot2(x,y)]
  20               [initials = .1:.2]
  21               [time= 0.0:40.0] [step = 40.0/200.0]
  22            };
  23           }; plot sol;
  24
  25
  26
  27
  28 # plot solution of
  29 #         [.1  .1  .2]
  30 #   X' =  [.1  2  .3] X  with initial condition t = 1, X = [1,2,3]
  31 #         [.22  1  1]
  32 function definitions {
  33           xdot1(x,y,z)= .1*x  +  .1*y  + .2*z;
  34           xdot2(x,y,z)= .1*x  + 2*y  +   .3*z;
  35           xdot3(x,y,z)= .22*x+   y  +   z;
  36          };
  37    sol = object{
  38         ode { [xdot1(x,y,z),xdot2(x,y,z),xdot3(x,y,z)]
  39               [initials = 1.0:2.0,3.0]
  40               [time= 1.0:5.0] [step = 4.0/200.0]
  41            };
  42 # plot (x,y,z) parametrically as t varies from 1.0 to 5.0
  43           }; plot sol;
  44
  45
  46
  47      A 2D Cantor set
  48
  49 set title "Triangular Cantor Set"
  50
  51         quatroot3=sqrt(3)*0.25;
  52         g11(x,y) = 0.5*x;     g12(x,y)=0.5*y;
  53         g21(x,y)=0.5*x+0.5;   g22(x,y)=0.5*y;
  54         g31(x,y)=0.5*x+0.25;  g32(x,y)=0.5*y+quatroot3;
  55
  56 Gx(x,y,z)=(z<=0.33? g31(x,y):(z<=0.67? g11(x,y):g21(x,y)));
  57 Gy(x,y,z)=(z<=0.33? g32(x,y):(z<=0.67? g12(x,y):g22(x,y)));
  58 Gz(x,y,z)=0.01*(rand()%100);
  59
  60 triangle = map{[Gx(x,y,z),Gy(x,y,z),Gz(x,y,z)]
  61                 [initial=0.:0.:1][itera=10000][color=green]}; plot2d triangle
  62
  63
  64     Solution to Burgers equation:
  65
  66 set dummy x,t;
  67 c1 = 0.5; c2 = 1.0;
  68
  69 f(x,t) = (c1*exp(c1*(x+c1*t)) +c2*exp(c2*(x+c2*t)))/ \
  70         (1. + exp(c1*(x+c1*t)) + exp(c2*(x+c2*t)));
  71
  72 a = sur{[8.0*f(x,t)][x=-10:20][t=-20:10]};
  73
  74 cup:
  75
  76
  77 r(x) = (x <= 0.0? -4*x + 0.1 : (x<= 2.0? 0.1: sqrt(  (1.3*(x-2)))))
  78
  79 glass = tube{[0,0,x][x=-0.3:4.5][radius=r(x)][sample=30:30]}
  80
  81
  82 torus:
  83
  84
  85 orusX(x) = cos(4*x) * (10.0 + 7.0 * cos(3*x))
  86 orusY(x) = sin(4*x) * (10.0 + 7.0 * cos(3*x))
  87 orusZ(x) = 7.0*sin(3*x)
  88
  89 tbe = tube{[orusX(x),orusY(x),orusZ(x)]
  90         [x= -pi:pi][samp=500:12]};
  91
  92
  93
  94 kdv:
  95 #-----------------------------------------------------
  96 # Simple Soliton soln to the KdV eq.
  97 #
  98
  99
 100 set title "Soliton Solution of the KdV Equation";
 101 set dummy u,v;
 102 c= 3.0/5.0;
 103 f1(u,v) = exp(c*u - c^3 * v);
 104 f2(u,v) = 1.0 + f1(u,v);
 105 kdv(u,v) = 3*c^2 * f1(u,v) / (f2(u,v))^2;
 106
 107 soliton = sur{[15*kdv(u,v)][u=-10:10][v=-10:10]};
 108
 109
 110 kdv2:
 111 #
 112 # Two Soliton soln to the KdV eq.
 113 #
 114
 115
 116 set dummy x,t;
 117 c1 = 0.5; c2 = 4./5.;
 118 e1(x,t)=exp(c1*x - c1^3 *t);
 119 e2(x,t)=exp(c2*x - c2^3 *t);
 120 a12(x,t) = ( (e1(x,t) - e2(x,t))/(e1(x,t)+e2(x,t)))^2;
 121 e12(x,t) = a12(x,t)*e1(x,t)*e2(x,t);
 122 f(x,t) = 1. + e1(x,t) + e2(x,t) + e12(x,t);
 123 fx(x,t) = c1 *e1(x,t) + c2 *e2(x,t) + (c1+c2)*e12(x,t);
 124 fc(x,t) = c1*e1(x,t) + c1 * e12(x,t);
 125 fxc(x,t) = c1^2*e1(x,t) + c1*(c1+c2)*e12(x,t);
 126 func(x,t) = (fxc(x,t)*f(x,t) - fx(x,t) *fc(x,t))/(f(x,t)^2);
 127
 128 a=sur{[30.*func(x,t)][x=-15:15][t=-15:15][sample=60:60]};
 129
 130
 131 map-hami:
 132 #
 133 # Plot 10 orbits of the hamiltonian map defined be
 134 #     [x,y]-->[hamX(x,y),hamY(x,y)]
 135 #---------------------------------------------------------
 136
 137
 138 function list {
 139         alpha = 0.211434654
 140         phi = 2.0*pi*alpha
 141         hamX(x,y) = x*cos(phi)-(y - x* x)*sin(phi)
 142         hamY(x,y) = x*sin(phi)+(y - x* x)*cos(phi)
 143         };
 144
 145 aa = object {
 146         map{[hamX(x,y),hamY(x,y)] [initial=0.16:0.578][itera=1000][color=4]}
 147         map{[hamX(x,y),hamY(x,y)] [initial=0.26:0.578][itera=1000][color=5]}
 148         map{[hamX(x,y),hamY(x,y)] [initial=0.268:0.578][itera=1000][color=10]}
 149         map{[hamX(x,y),hamY(x,y)] [initial=0.38:0.578][itera=1000][color=9]}
 150         map{[hamX(x,y),hamY(x,y)] [initial=0.36:0.478][itera=1000][color=6]}
 151         map{[hamX(x,y),hamY(x,y)] [initial=0.46:0.378][itera=1000][color=3]}
 152         map{[hamX(x,y),hamY(x,y)] [initial=0.06:0.578][itera=1000][color=2]}
 153         map{[hamX(x,y),hamY(x,y)] [initial=0.46:0.300][itera=1000][color=1]}
 154         map{[hamX(x,y),hamY(x,y)] [initial=0.20:0.278][itera=1000][color=7]}
 155         map{[hamX(x,y),hamY(x,y)] [initial=0.10:0.178][itera=1000][color=8]}
 156         }; plot2d aa
 157
 158
 159 map-quad:
 160
 161
 162 lambda = 3.5;
 163 f(x) = lambda * x*(1-x);
 164
 165 a = cur{[x,f(x),0][x=0:1][samples=300][color=red]};  # y = f(x)
 166 b = cur{[x,x,0] [x=0,1][sample=10][color=green]};    # y = x
 167
 168 # the map
 169 #      (x, y) -> (y, f(x))
 170 # The z-component is not relevent (in plot2d) but it is
 171 # helpful to get a depthcue shading.
 172 #
 173
 174 c = map{[y,f(x),0.002+z][initial=0.1,f(0.1),-1]
 175         [iterates=400][color=blue][style=solid]}
 176
 177 plot2d  a, b, c;
 178
 179
 180 minimal:
 181
 182 set dummy u,v,alpha, beta;
 183 #
 184 #  Enneper's surface
 185 #
 186 Enneper = sur{[u-u**3/3+u*v**2,v-v**3/3+v*u**2,u**2-v**2]
 187                 [u=-1:1][v=-1:1]}
 188
 189
 190 minimal1:
 191
 192 Catanoid = sur{[cosh(u)*cos(v), 3*u, -cosh(u)*sin(v)]
 193                 [u=-pi:pi][v=-pi:pi]};
 194
 195
 196
 197 minimal2:
 198
 199 set dummy u, v;
 200
 201 X(u,v) = cosh(2*u)*cos(2*v) -1.0;
 202 Y(u,v) = -sinh(u)*sin(v) - 1./3. *sinh(3*u)*sin(3*v);
 203 Z(u,v) = -sinh(u)*cos(v) - 1./3. *sinh(3*u)*cos(3*v);
 204
 205 e1 =sur{[1.5*X(u,v),Y(u,v),Z(u,v)]
 206         [u = -1: 1][ v= -1: 1]
 207         [sample = 40, 40]};
 208
 209
 210 minimal1:
 211
 212 # Higher order Enneper's surfaces
 213 #
 214 set dummy u,v;
 215 set title "Enneper's minimal surface";
 216 n = 2;
 217
 218 if(n > 1) {set title "Higher order Enneper's minimal surface";}
 219
 220
 221 X(z) = real( 0.5*(z - z^(2*n+1)/(2*n+1)) );
 222 Y(z) = imag( 0.5*(z + z^(2*n+1)/(2*n+1)) );
 223 Z(z) = real( z^(n+1)/(n+1));
 224 e1 =sur{[X( (v)*(sin(u) + i* cos(u))),
 225          Y( (v)*(sin(u) + i *cos(u))),
 226          Z( (v)*(sin(u) + i* cos(u)))]
 227         [u = 0: pi][ v= -1.2: 1.2]
 228         [sample = 40*n, 40]};
 229
 230
 231 minimal3:
 232
 233 #
 234 set dummy u,v;
 235 n = 1;
 236
 237 X(z) = real( 0.5*(1./z - z^(2*n+1)/(2*n+1)) );
 238 Y(z) = imag( 0.5*(-1./z + z^(2*n+1)/(2*n+1)) );
 239 Z(z) = real( z^n/n);
 240
 241 e1 =sur{[X( v*sin(u) + i* v *cos(u)),
 242          Y( v*sin(u) + i* v *cos(u)),
 243          Z( v*sin(u) + i* v *cos(u))]
 244         [u = -pi: pi][ v= 1.0/3.0: 6.0/5.0]
 245         [sample = 80*n, 80]};
 246
 247
 248
 249 misc1:
 250
 251 set dummy u,v;
 252 set title "Yet Another Parametric Surface"
 253
 254 a=sur{[sin(u)*cos(v), sin(u)*sin(v),-0.25*cosh(u)]
 255         [u=0:pi][v=0:2*pi]};
 256
 257
 258 misc3:
 259
 260 set dummy u,v;
 261 set title "Yet Another Parametric Surface"
 262 a=sur{[u*cos(v)*sin(u), u*cos(u)*cos(v), -u*sin(v)]
 263         [u=0:3*pi][v=0:pi]};
 264
 265
 266 objs:
 267
 268 #   another example of multiple objects
 269 #
 270
 271 f(x) = sgn(sin(x))* ( (abs(sin(x)))^0.4);
 272 g(x) =  sgn(cos(x))* ( (abs(cos(x)))^0.4);
 273
 274 a=tube{[ f(x), g(x), -1.3]
 275         [x=0:2*pi][radius=0.1]
 276         [sample=40:6]}
 277 b=tube{[ f(x), g(x), 1.3]
 278         [x=0:2*pi][radius=0.1]
 279         [sample=40:6]}
 280
 281 c=tube{[ f(x), 1.3, g(x)]
 282         [x=0:2*pi][radius=0.1]
 283         [sample=40:6]}
 284 d=tube{[ f(x), -1.3, g(x)]
 285         [x=0:2*pi][radius=0.1]
 286         [sample=40:6]}
 287
 288 e=tube{[ -1.3 f(x), g(x)]
 289         [x=0:2*pi][radius=0.1]
 290         [sample=40:6]}
 291 f=tube{[ 1.3 f(x), g(x)]
 292         [x=0:2*pi][radius=0.1]
 293         [sample=40:6]}
 294
 295 o1 = sur{[cos(x)*cos(y),cos(x)*sin(y),sin(x)][x=-0.5*pi:0.5*pi]
 296          [y=0:2*pi][samples= 16:20]};
 297 o2=obj{a,b,c,d,e,f};
 298
 299 plot o1,o2;
 300
 301
 302 ode-lorenz:
 303
 304 #
 305 # Note: Need about 3 minutes to solve this equation
 306 #
 307 #---------------------------------------------------------
 308
 309 function definitions {
 310         RR = 28.0
 311         LORdy1dt (x,y,z) = 10.0* (y - x)
 312         LORdy2dt (x,y,z) = RR* x - x*z - y
 313         LORdy3dt (x,y,z) = x* y - 8.0* z /3.0
 314         };
 315
 316 lorenz = object{
 317              eqn{
 318                  [LORdy1dt(x,y,z),LORdy2dt(x,y,z),LORdy3dt(x,y,z)]
 319                  [initials = 0.03: 0.12: 0.1]
 320                  [time = 0.0:80.0 ][step = 0.3]
 321                 };
 322               }; plot lorenz;
 323
 324
 325
 326 ode-rossler:
 327
 328 # Note: Need about 10 minutes to solve this equation
 329 #
 330 #---------------------------------------------------------
 331
 332 ROS = 5.6
 333 ROSdy1dt (x,y,z) = -(y + z)
 334 ROSdy2dt (x,y,z) = x + 0.2 * y
 335 ROSdy3dt (x,y,z) = 0.2 + z * (x - ROS)
 336
 337 rossler = eqn{
 338       [ROSdy1dt(x,y,z),ROSdy2dt(x,y,z),ROSdy3dt(x,y,z)]
 339       [initials = 0.46:0.20:0.69]
 340       [time = 0.0:500.0]  [step = 0.002] [skip 5]
 341       [method = RKQC][small = 0.00000001]
 342       [axes: automatic]
 343     }; plot rossler;
 344
 345
 346
 347 star2:
 348
 349 a=9; b=6; c=6;
 350 f(x,y) = (2+sin(a*x)*0.5)*cos(x+sin(b*x)/c)*y;
 351 g(x,y) = (2+sin(a*x)*0.5)*sin(x+sin(b*x)/c)*y;
 352
 353 p2 = sur{[f(x,y),g(x,y),2*sin(x)*y]
 354          [x=0, 2*pi][y=0:1][sample= 300:10]};
 355
 356 plot2d p2
 357
 358
 359 sur-sqrtz:
 360
 361 # The Riemann surface for w = sqrt(z) is realized by picing the
 362 # 2 branches together
 363
 364 zz(x,y) = x *exp(i * y)
 365 sqtzX(x,y) = real(zz(x,y))
 366 sqtzY(x,y) = imag(zz(x,y))
 367 sqtzZ(x,y) = real( sqrt(zz(x,y)))
 368
 369 # we need to define 2 branches, one with sqrt(1) = 1; the other with
 370 # sqrt(1) = -1
 371
 372 sqrt = object{
 373         mat default;
 374         sur{[sqtzX(x,y),sqtzY(x,y),sqtzZ(x,y)] [x=0:1][y=0:2*pi] [samp=32:32]}
 375         sur{[sqtzX(x,y),sqtzY(x,y),-sqtzZ(x,y)] [x=0:1][y=0:2*pi] [samp=32:32]}
 376         }; plot sqrt
 377
 378
 379 sur-mobius:
 380
 381 #                ----------- The mobuis band.
 382 #
 383 #---------------------------------------------------------------
 384
 385 MobuisX(x,y) = cos(x) * (2.0 + y * cos(0.5 * x))
 386 MobuisY(x,y) = sin(x) * (2.0 + y * cos(0.5 * x))
 387 MobuisZ(x,y) = y * sin(0.5 * x)
 388
 389 mobius = surf{
 390               [MobuisX(x,y),MobuisY(x,y),MobuisZ(x,y)]
 391               [x=-pi:pi][y=-0.5:0.5][samp=60:10]
 392             };
 393 #---------------------------------------------------------------
 394
 395
 396 sur-klein:
 397
 398 set dummy u,v;
 399 set title "The Klein Bottle";
 400
 401
 402 KleinX(u,v) = (2*sin(u)*cos(v/2)-sin(2*u)*sin(v/2)+8)*cos(v);
 403 KleinY(u,v) = (2*sin(u)*cos(v/2)-sin(2*u)*sin(v/2)+8)*sin(v);
 404 KleinZ(u,v) = 2*sin(u)*sin(v/2)+sin(2*u)*cos(v/2);
 405
 406 klein = surf{
 407            [KleinX(u,v),KleinY(u,v),KleinZ(u,v)]
 408            [u=-pi:pi][v=-pi:pi][samp=50:30]
 409           };
 410
 411
 412
 413 tube:
 414
 415 set title " A Torus Knot"
 416
 417 rusX(x) = cos(3*x) * (10.0 + 6.0 * cos(2*x))
 418 rusY(x) = sin(3*x) * (10.0 + 6.0 * cos(2*x))
 419 rusZ(x) = -6.0*sin(2*x)
 420
 421 a = tube{[rusX(x),rusY(x),rusZ(x)]
 422       [x= -pi:pi][samp=200:20][radius=1.4 + 0.1*sin(16*x)]};
 423
 424 a = tube{[rusX(x),rusY(x),rusZ(x)]
 425       [x= -pi:pi][samp=200:20][radius=1.4][sample=100:16]};
 426
 427
 428
 429 xplot:
 430
 431 graph list{
 432         x1 = tube{[15+x,1.7*x,-5][x=-7:7][sample=30:10]
 433                 [radius=0.6*cos(3*x)+1.2]};
 434         x2 = tube{[15+x, -1.7*x,-5][x=-7:7][sample=30:10]
 435                 [radius=0.6*cos(3*x)+1.2]};
 436         p1 = tube{[5.8, x,0][x=-4:4][sample=3:10][radius=0.7]};
 437         p2 = tube{[6.4+2.8*cos(x),1.5-2*sin(x),0][x=-0.5*pi,0.5*pi]\
 438                 [sample=10:10][radius=0.7]};
 439         l1 = tube{[11,x,0][x=-4:4][sample=2:10][radius=0.7]};
 440         l2 = tube{[11+x,-3.5,0][x=0,3.5][sample=2:10][radius=0.7]};
 441         o =  tube{[17.5+1.8*sgn(cos(x))* ( abs(cos(x)))^0.4,
 442                         3*sgn(sin(x))* ( abs(sin(x)))^0.4,0]
 443                 [x=0:2*pi][sample=20:10] [radius=0.7]};
 444         t1 = tube{[20.8+x,3.0,0][x=0:5][sample=2:10][radius=0.7]};
 445         t2 = tube{[23.3,x,0][x=-4:2.5][sample=2:10][radius=0.7]};
 446
 447 };
 448 X = obj{x1,x2};
 449 P = obj{p1,p2};
 450 L = obj{l1,l2};
 451 O = obj{o};
 452 T = obj{t1,t2};
 453
 454 plot2d  X P L O T;
 455
 456
 457