archive/books/schelter/427k-2.bk

   1 \x06\x01\x19\x16\x05
   2 ((face (dfplot-eval 3029 3093 3544 3634) (maxima-eval 1513 1528 2348 2398 2854 2921) (book-result 554 560 589 590 630 675 757 824 870 922 976 1068 1619 1645 1732 1759 2191 2281 2459 2549) (maxima-eval-insert 531 545 566 575 597 607 703 746 831 860 926 965 1690 1723 2105 2181 2405 2449) (bold 7 93 1763 1778)) (book-command-arg))\f
   3
   4     Tutorial on solving some 1st and 2nd order ODES using Symbolic
   5                   Computation.
   6
   7
   8 In this section we use maxima to obtain symbolic solutions to certain
   9 differential equations.  Then at the end we plot these, and
  10 also compare with the dfplot trajectories.
  11
  12 We use the ode2 routine:
  13
  14 -- Function ode2( differential_equation, dependent_var, independent_var )
  15
  16 To enter a differential, you must put a quote "'" before the diff.
  17 This quote stops the derivative from actually being evaluated.
  18
  19   diff(sin(x),x) returns COS(X) and  diff(y,x) evaluates to 0
  20
  21 but 'diff(y,x) yields what we want:
  22
  23                                       dY
  24                                       --
  25                                       dX
  26
  27 [see notation below].
  28
  29     ode2(x^2*'diff(y,x)+3*x*y = sin(x)/x,  y,x) returns
  30
  31
  32                                     %C - COS(X)
  33                                 Y = -----------
  34                                          3
  35                                         X
  36
  37
  38
  39    ode2('diff(y,x)+y^3 = 0 ,y,x) yields
  40
  41
  42                                   1
  43                                  ---- = X + %C
  44                                     2
  45                                  2 Y
  46
  47
  48   ode2('diff(y,x) = 2*(1+x)*(1+y^2) ,y,x) returns
  49
  50
  51                                        2
  52                             ATAN(Y)   X  + 2 X
  53                             ------- = -------- + %C
  54                                2         2
  55
  56
  57 Where is the above solution valid?
  58
  59 Notation: In maxima variables may have multiple letters, so that 'ab'
  60 is one variable, and you must write 'a*b' to mean the product of two
  61 symbols named 'a' and 'b'.  Similarly 'dy' is just a two letter
  62 symbol.  The printing of dy/dx is a pretty printing for 'diff(y,x) ,
  63 but you must use this latter representation to enter the form.  To get
  64 a display which uses the same notation for input and output do
  65 display2d:false .  If you do that the last form evals to the more
  66 compact but perhaps less human readable ATAN(Y)/2 = (X^2+2*X)/2+%C .
  67
  68 A second order differential equation:
  69
  70   ode2('diff(y,x,2) = 2*(1+x) ,y,x) returns Y = (X^3+3*X^2)/3+%K2*X+%K1
  71
  72
  73 Initial Values:
  74
  75 It is possible to have the solver substitute in initial values, to determine
  76 the constant(s), in the above solutions.   The functions for doing that
  77 are
  78
  79 IC1(sol,x=x0,y=y0)
  80
  81 IC2 is used for specifying an initial tangent and a point to pass through:
  82 For example the following forces the slope to be 3 at the point (1,2).
  83
  84     (soln:ode2('diff(y,x,2) = 2*(1+x) ,y,x), tmp:IC2(soln,x=1,y=2,'diff(y,x)=3)) returns
  85
  86                                     3      2
  87                                    X  + 3 X    2
  88                                Y = --------- + -
  89                                        3       3
  90
  91
  92   To plot this, we had set tmp to be the solution equation,
  93
  94    plot2d(sublis(solve(tmp,y),y),[x,-4,4],[y,-10,10])
  95
  96
  97
  98   sol:ode2('diff(y,x) = 2*(1+x)*(1+y**2) ,y,x) returns
  99
 100                                        2
 101                             ATAN(Y)   X  + 2 X
 102                             ------- = -------- + %C
 103                                2         2
 104
 105    Lets plot this solution for x in [2,4] and y in [-3,6], and verify
 106 that the point (3,4) does lie on the particular solution computed by
 107 IC1.  After you have clicked on the above, to set sol correctly we can
 108 plot the particular solution (in order to keep plot2d happy, we have
 109 to solve for y) :
 110
 111       plot2d(  sublis(solve( IC1(sol,x=3,y=4) ,y),y),   [x,2,4],[y,-3,6])
 112
 113  We could look at the following, and click at the position (3,4) to
 114 find the solution numerically.
 115
 116      ode{d[y,x]=2*(1+x)*(1+y**2)};set xrange [2,4];set yrange[-10,10]
 117
 118 This uses the default 4th order Runge Kutte with a fixed step in the
 119 independent variable x.  Note that this gives an answer which is
 120 sloped a bit too much to the right at the top, whereas our theoretical
 121 plot was practically vertical when going through (3,4).  We can retry
 122 using an 'adaptive 5th order Runge Kutta method (RKQC)'.  This alters
 123 the x step according to an error estimate, and gives a more accurate
 124 picture.  It is also slower.
 125
 126      ode{d[y,x]=2*(1+x)*(1+y**2), [method=RKQC]};
 127           set xrange [2,4];set yrange[-10,10]
 128
 129
 130
 131
 132
 133
 134
 135
 136
 137