macro(h_cc = 1); macro(m_cc = 2); macro(its_cc = 3); macro(steps_cc = 4); macro(ti_cc = 5); macro(tf_cc = 6); macro(intmethod_cc = 7); macro(Reduct_cc = 8); macro(a_cc = 9); macro(res_cc = 10); macro(Idid_cc = 11); macro(numdiff_cc = 12); macro(method_cc = 13); macro(sUse_cc = 14); macro(ss_cc = 15); macro(hmin_cc = 16); macro(hmax_cc = 17); macro(Direc_cc = 18); macro(maxIt_cc = 19); macro(Ecur_cc = 20); macro(polord_cc = 21); macro(laststored_cc = 22); macro(AbsTol_cc = 23); macro(RelTol_cc = 24); macro(tout_cc = 25); macro(sUsed_cc = 26); macro(Emax_cc = 27); macro(sMax_cc = 28); macro(sOld_cc = 29); macro(sStep_cc = 30); macro(p_cc = 31); macro(q_cc = 32); macro(r_cc = 33); macro(s_cc = 34); macro(hboxes_cc = 35); macro(vboxes_cc = 36); macro(nclines_cc = 37); macro(flowf_cc = 38); macro(segments_cc = 39); macro(orbthick_cc = 40); macro(asp_cc = 41); macro(fhead_cc = 42); macro(fishs_cc = 43); macro(fishl_cc = 44); macro(u_cc = 45); macro(v_cc = 46); macro(fthick_cc = 47); macro(field_cc = 48); macro(rkf45_cc = 4); macro(gear_cc = 5); macro(mgear_cc = 6); macro(gear_cc = 7); macro(nreach = 1); macro(fail = 4); macro(singul = 5); macro(userv = 6); macro(nexth = 7); macro(ddone = 8); macro(polerr = 9); macro(node = 10); macro(_AA = `plots/ODE/AA`); macro(_body = `plots/ODE/body`); macro(_eq = `plots/ODE/eq`); macro(gray = COLOR(RGB, 0.5,0.5,0.5)): macro(black = COLOR(RGB, 0,0,0)): macro(red = COLOR(RGB, 1, 0, 0)); `besirk/sSeq` := array(1..13, [1, 2, 3, 5, 7, 11, 13, 17, 23, 31, 47, 61, 89]): alias(redscale = (xx -> COLOR(RGB,1,1-xx,1-xx))): alias(rungekuttahf = rungekutta): `ODE/setup` := proc(F, CC, optargs) local inits, carg, initlist, vars, precision, m, i, tt, FF, ttt, oeq, oneq, func, j, vv, intm, indep; precision := max(Digits, trunc(evalhf(Digits))); intm := subs(select(xx -> (op(1, xx) = 'intmethod'),optargs), 'intmethod'); if intm = 'intmethod' then CC[intmethod_cc] := 0 elif intm = 'besirk' then CC[intmethod_cc] := 1 elif intm = 'dsolve' then CC[intmethod_cc] := 2 else CC[intmethod_cc] := 0 fi; initlist := select(xx -> type(rhs(xx),{set,list}), optargs); initlist := select(xx -> (lhs(xx) <> 'grid' and lhs(xx) <> 'view' and lhs(xx) <> 'variables' and lhs(xx) <> 'plotrange' and lhs(xx) <> 'projections'), initlist); if nops(initlist) <> 0 then inits := rhs(op(initlist)); if type(inits,set) then inits := [op(inits)] elif type(inits[1],list) then inits := inits else inits := [inits] fi; if CC[m_cc] = 2 then if nops(inits[1]) = 2 then inits := map( x -> [0,op(x)], inits) fi fi; CC[m_cc] := nops(inits[1])-1 else inits := [] fi; m := CC[m_cc]; if type(F, procedure) then vars := [op(1,eval(F))]; if m = 1 then func := array([F(_AA[1],_AA[2])]) else func := array(F(seq(_AA[i],i=1..m+1))); fi; elif type(F[1], procedure) then func := array([seq(F[j](seq(_AA[i],i=1..m+1)),j=1..m)]); vars := [op(1,eval(F[1]))]; else if not type(F,list) then FF := [F] else FF := F fi; vars := subs(optargs, variables); if vars = 'variables' then vars := [op(2,lhs(FF[1]))]; for i to m do vars := [op(vars), op(1,lhs(FF[i]))]; od; fi; func := array(subs({seq(vars[i] = _AA[i] , i=1..m+1)}, [seq(rhs(FF[i]),i=1..m)])); vars := [vars[1], op(map(proc(unk, ut) if evalb(op(1, unk) = ut) then op(0, unk) else unk fi end, [seq(vars[i],i=2..m+1)], vars[1]))]; fi; tt := select(xx -> type(rhs(xx), range), optargs); vv := subs(select(xx -> (op(1, xx) = 'view'),optargs), 'view'); if nops(inits) > 0 then ttt := subs({op(tt)}, vars[1]); if not type(ttt, range) then ttt := subs({op(tt)}, 'tdummy') fi; if not type(ttt, range) then if not vv = 'view' then ttt := subs(vv, vars[1]) fi fi; if not type(ttt, range) then if nops(tt) = 0 then ttt := -2.5 .. 2.5 else if m = 1 then tt := subsop(1 = (vars[1] = rhs(tt[1])), tt) else tt := subsop(nops(tt) = (vars[1] = rhs(tt[nops(tt)])), tt) fi; ttt := subs({op(tt)}, vars[1]) fi fi; CC[ti_cc] := op(1, ttt); CC[tf_cc] := op(2, ttt) fi; if nops(initlist) <> 0 then if round(CC[intmethod_cc]) = 0 or round(CC[intmethod_cc]) = 2 then `ODE/rkhf/setcc`(F, CC, optargs); elif round(CC[intmethod_cc]) = 1 then `besirk/setcc`(F, CC, optargs) fi fi; [inits, vars, func] end: `ODE/graphics/setcc` := proc(F, CC, optargs, vars) local carg, precision, tt, vv, m, vvv, k, i, ttt; m := CC[m_cc]; CC[nclines_cc] := 0; CC[flowf_cc] := 0; precision := max(Digits, trunc(evalhf(Digits))); CC[segments_cc] := 3; CC[asp_cc] := 1; CC[hboxes_cc] := 14: CC[vboxes_cc] := 14: CC[fhead_cc] := 1; CC[fishs_cc] := 1; CC[fishl_cc] := 1; CC[orbthick_cc] := 1; CC[fthick_cc] := 1; for carg in optargs do if not type(carg, equation) then ERROR(`optional arguments must be equations`) elif op(1, carg) = 'flowthickness' then CC[fhead_cc] := op(2, carg); elif op(1, carg) = 'fieldthickness' then CC[fthick_cc] := op(2, carg); elif op(1, carg) = 'orbitthickness' then CC[orbthick_cc] := op(2, carg); elif op(1, carg) = 'flowlength' then CC[fishs_cc] := op(2, carg); elif op(1, carg) = 'flowlogscale' then CC[fishl_cc] := op(2, carg); elif op(1, carg) = 'segments' then CC[segments_cc] := max(3,op(2, carg)); elif op(1, carg) = 'aspectratio' then CC[asp_cc] := op(2, carg); elif op(1, carg) = 'grid' then if not type(rhs(carg), [integer, integer]) then ERROR(`grid must be a list of two positive integers`) fi; CC[hboxes_cc] := op(2, carg)[1]; CC[vboxes_cc] := op(2, carg)[2]; elif op(1, carg) = 'parametricplot' then if op(2, carg) = 'flowparametricplot' then CC[flowf_cc] := 1 elif op(2, carg) = 'vectorfieldplot' then CC[flowf_cc] := 2 else CC[flowf_cc] := 0 fi elif op(1, carg) = 'directionfield' then if op(2, carg) then CC[field_cc] := 1 else CC[field_cc] := 0 fi elif op(1, carg) = 'vectorfield' then if op(2, carg) then CC[field_cc] := 2 else CC[field_cc] := 0 fi elif op(1, carg) = 'flowfield' then if op(2, carg) then CC[field_cc] := 3 else CC[field_cc] := 0 fi elif op(1, carg) = `nullclines` then if op(2, carg) then CC[nclines_cc] := 1 fi; elif op(1, carg) = `isoclines` then if op(2, carg) then CC[nclines_cc] := 1 fi fi od; vv := subs(select(xx -> (op(1, xx) = 'view'),optargs), 'view'); tt := select(xx -> type(rhs(xx), range), optargs); if vv = 'view' then vvv := subs({op(tt)}, [seq(vars[k], k = 1 .. m + 1)]); vv := []; for i to m + 1 do vv := [op(vv), vars[i] = vvv[i]] od; vv := select(xx -> type(rhs(xx), range), vv); if nops(vv) < 2 then vv := [seq([vars[1],vars[k]], k = 2..m+1)] else vv := [map(xx -> lhs(xx), vv)] fi elif not type(vv[1],list) then vv := [vv] fi; ttt := subs({op(tt)}, [seq(vv[1][i], i = 1 .. nops(vv[1]))]); if type(ttt, {[range, range], [range, range, range]}) then CC[p_cc] := op(1, ttt[1]); CC[q_cc] := op(2, ttt[1]); CC[r_cc] := op(1, ttt[2]); CC[s_cc] := op(2, ttt[2]); CC[u_cc] := op(1, ttt[nops(ttt)]); CC[v_cc] := op(2, ttt[nops(ttt)]) else ttt := [seq('DEFAULT', i = 1 .. nops(vv[1]))]; fi; [vv, ttt] end: `ODE/rkhf/setcc` := proc(F, CC, optargs) local tt, initlist, m, precision, carg; precision := max(Digits, trunc(evalhf(Digits))); CC[hmin_cc] := 10^(-precision+3); CC[steps_cc] := 50; CC[res_cc] := 0; CC[its_cc] := 1; CC[h_cc] := evalf((CC[tf_cc] - CC[ti_cc])/50, precision); for carg in optargs do if not type(carg, equation) then ERROR(`optional arguments must be equations`) elif op(1, carg) = 'stepsize' then CC[h_cc] := evalf(op(2, carg), precision); CC[steps_cc] := round(evalf((CC[tf_cc] - CC[ti_cc])/op(2, carg), precision)); elif op(1, carg) = 'bound' then if op(2, carg) then CC[res_cc] := 1 fi; elif op(1, carg) = 'iterations' then CC[its_cc] := evalf(op(2, carg), precision); elif op(1, carg) = 'numsteps' then CC[steps_cc] := op(2,carg); CC[h_cc] := evalf((CC[tf_cc] - CC[ti_cc])/op(2, carg), precision) elif op(1, carg) = 'method' then if op(2, carg) = 'rkf45' then CC[method_cc] := rkf45_cc elif op(2, carg) = 'gear' then CC[method_cc] := gear_cc elif op(2, carg) = 'mgear' then CC[method_cc] := mgear_cc elif op(2, carg) = 'dverk78' then CC[method_cc] := gear_cc fi fi; od; end: `besirk/setcc` := proc(F, CC, optargs) local m, precision, carg; precision := max(Digits, trunc(evalhf(Digits))); m := CC[m_cc]; CC[numdiff_cc] := `ODE/getnumdiff`(F, optargs, m); CC[hmax_cc] := evalf((CC[tf_cc] - CC[ti_cc])/50, precision); CC[h_cc] := CC[hmax_cc]; CC[a_cc] := evalf(-1/2*2^(1/2)*cos(1/3*arctan(1/4*2^(1/2))) + 1 + 1/2* 3^(1/2)*2^(1/2)*sin(1/3*arctan(1/4*2^(1/2))), precision); CC[method_cc] := 1; CC[maxIt_cc] := 1000; CC[Emax_cc] := 2; CC[sMax_cc] := 12; CC[res_cc] := 0; CC[hmin_cc] := .0001; CC[sUse_cc] := 11; CC[polord_cc] := 3; CC[AbsTol_cc] := 10^(7 - precision); CC[RelTol_cc] := 0; CC[Idid_cc] := 0; for carg in optargs do if not type(carg, equation) then ERROR(`optional arguments must be equations`) elif op(1, carg) = 'numsteps' then CC[steps_cc] := op(2,carg); CC[hmax_cc] := evalf((CC[tf_cc] - CC[ti_cc])/op(2, carg), precision); CC[h_cc] := CC[hmax_cc] elif op(1, carg) = 'hmin' then CC[hmin_cc] := evalf(op(2, carg), precision) elif op(1, carg) = 'hmax' then CC[hmax_cc] := evalf(op(2, carg), precision); CC[h_cc] := CC[hmax_cc] elif op(1, carg) = 'stepsize' then CC[hmax_cc] := evalf(op(2, carg), precision); CC[h_cc] := CC[hmax_cc] elif op(1, carg) = 'maxsteps' then CC[maxIt_cc] := op(2, carg) elif op(1, carg) = 'AbsTolerance' then CC[AbsTol_cc] := op(2, carg) elif op(1, carg) = 'RelTolerance' then CC[RelTol_cc] := op(2, carg) elif op(1, carg) = 'polord' then CC[polord_cc] := op(2, carg) elif op(1, carg) = 'bound' then if op(2, carg) then CC[res_cc] := 1 fi; elif op(1, carg) = 'startstepsize' then CC[h_cc] := evalf(op(2, carg), precision) elif op(1, carg) = 'method' then if op(2, carg) = 'rk4' then CC[method_cc] := 2 elif op(2, carg) = 'rk2' then CC[method_cc] := 3 fi fi od; end: firsteuler := proc(f, initlist, h0, n, iterations) local i, xk, yk, pts, h, npts; options `Copyright 1991 by Daniel Schwalbe`; pts := array(0..n); if nargs > 4 then npts := iterations; else npts := 1 fi; h := evalf(h0/npts); xk := evalf(initlist[1]); yk := evalf(initlist[2]); pts[0] := [xk, yk]; for i to n while type(yk, numeric) do to npts while type(yk, numeric) do yk := yk + h*f(xk, yk); xk := xk + h od; pts[i] := [xk, yk] od: if i <= n then ERROR(`undefined procedure`) else op(pts) fi: end: impeuler := proc(f, initlist, h0, n, iterations) local h, xk, yk, pts, j, npts: options `Copyright 1991 by Daniel Schwalbe`; pts := array(0..n); if nargs > 4 then npts := iterations else npts := 1 fi; h := evalf(h0/npts): xk := evalf(initlist[1]); yk := evalf(initlist[2]); pts[0] := [xk, yk]; for j to n while type(yk, numeric) do to npts while type(yk, numeric) do yk := yk + h/2*(f(xk, yk) + f(xk + h, yk + h*f(xk, yk))): xk := xk + h od: pts[j] := [xk, yk] od: if j <= n then ERROR(`undefined procedure`) else op(pts) fi: end: `ODE/rkhf` := proc(f, B, CC) local j, i, h, n, npts, t, k1, k2, k3, k4, m, y: options `Copyright 1993 by Daniel Schwalbe`; h := CC[h_cc]; n := CC[steps_cc]; npts := CC[its_cc]; m := CC[m_cc]; t := array(1..m + 1); y := array(1..m + 1); for i to m + 1 do t[i] := B[0, i - 1] od; k1 := array(1..m);k2 := array(1..m); k3 := array(1..m);k4 := array(1..m); for j from 0 to n - 1 do to npts do f(t, k1): y[1] := t[1] + h/2; for i from 2 to m + 1 do y[i] := t[i] + h*k1[i - 1]/2 od; f(y, k2): for i from 2 to m + 1 do y[i] := t[i] + h*k2[i - 1]/2 od; f(y, k3): y[1] := t[1] + h; for i from 2 to m + 1 do y[i] := t[i] + h*k3[i - 1] od; f(y, k4): t[1] := y[1]; for i from 2 to m + 1 do t[i] := t[i] + h*(k1[i - 1] + 2*k2[i - 1] + 2*k3[i - 1] + k4[i - 1])/6 od od: for i to m + 1 do B[j + 1, i - 1] := t[i] od: od: end: `ODE/rkhf/bounded` := proc(f, B, CC) local j, i, h, n, npts, t, k1, k2, k3, k4, m, y, v,p,q,r,s, dy, hm, hmin: options `Copyright 1993 by Daniel Schwalbe`; h := CC[h_cc]; n := CC[steps_cc]; npts := CC[its_cc]; m := CC[m_cc]; v := round(min(m + 1, 3)); p := CC[p_cc]; q := CC[q_cc]; r := CC[r_cc]; s := CC[s_cc]; t := array(1..m + 1); y := array(1..m + 1); for i to m + 1 do t[i] := B[0, i - 1] od; k1 := array(1..m);k2 := array(1..m); k3 := array(1..m);k4 := array(1..m); hm := CC[hmax_cc]; hmin := CC[hmin_cc]; for j from 0 to n - 1 do to npts do f(t, k1): y[1] := t[1] + h/2; for i from 2 to m + 1 do y[i] := t[i] + h*k1[i - 1]/2 od; f(y, k2): for i from 2 to m + 1 do y[i] := t[i] + h*k2[i - 1]/2 od; f(y, k3): y[1] := t[1] + h; for i from 2 to m + 1 do y[i] := t[i] + h*k3[i - 1] od; f(y, k4): t[1] := y[1]; for i from 2 to m + 1 do y[i] := h*(k1[i - 1] + 2*k2[i - 1] + 2*k3[i - 1] + k4[i - 1])/6; t[i] := t[i] + y[i] od; if abs(h) > 0 and ((t[v] < r or t[v] > s or t[v-1] < p or t[v-1] > q) or (abs(t[v] - B[j, v-1]) < hmin and abs(t[v-1] - B[j, v-2]) < 10^(-7))) then CC[steps_cc] := j+1; h := 0 else hm := max(hm, abs(y[2])) fi; od: for i to m + 1 do B[j + 1, i - 1] := t[i] od: od: CC[hmax_cc] := hm; end: `ODE/rkhfbandf` := proc(f, CC, init) local nf, nb, precision, m, Af, Ab, B, i, k, n, h; precision := max(Digits, trunc(evalhf(Digits))); m := round(CC[m_cc]); h := CC[h_cc]; nf := round(evalf((CC[tf_cc] - init[1])/h)); nb := round(evalf((init[1] - CC[ti_cc])/h)); Af := array(0..nf, 0..m ); Ab := array(0..nb, 0..m ); for i to m + 1 do Af[0, i - 1] := evalf(init[i], precision) od; for i to m + 1 do Ab[0, i - 1] := evalf(init[i], precision) od; n := CC[steps_cc]; CC[h_cc] := h/CC[its_cc]; CC[steps_cc] := nf; if round(CC[res_cc]) = 0 then if precision <= evalhf(Digits) then evalhf(`ODE/rkhf`(f, var(Af), var(CC))) else `ODE/rkhf`(f, Af, CC) fi; else if precision <= evalhf(Digits) then evalhf(`ODE/rkhf/bounded`(f, var(Af), var(CC))) else `ODE/rkhf/bounded`(f, Af, CC) fi; fi; nf := round(min(nf, (trunc(CC[steps_cc]/CC[segments_cc])+1)*CC[segments_cc])); CC[steps_cc] := nb; CC[h_cc] := -CC[h_cc]; if round(CC[res_cc]) = 0 then if precision <= evalhf(Digits) then evalhf(`ODE/rkhf`(f, var(Ab), var(CC))) else `ODE/rkhf`(f, Ab, CC) fi; else if precision <= evalhf(Digits) then evalhf(`ODE/rkhf/bounded`(f, var(Ab), var(CC))) else `ODE/rkhf/bounded`(f, Ab, CC) fi; fi; nb := round(min(nb, (trunc(CC[steps_cc]/CC[segments_cc])+1)*CC[segments_cc])); CC[h_cc] := h; CC[steps_cc] := n; B := array(-nb .. nf); for i from 0 to nf do B[i] := [seq(Af[i, k - 1], k = 1..m + 1)] od; for i to nb do B[-i] := [seq(Ab[i, k - 1], k = 1..m + 1)] od; eval(B) end: `ODE/dsolve` := proc(func, CC, init) local nf, nb, m, B, i, j, h, ff, t, yy, v, mt; m := round(CC[m_cc]); h := CC[h_cc]; B := [init]; v := seq(convert(v[i], string), i = 1..m+1); ff := subs({seq(_AA[i] = convert(v[i],string),i=1..m+1)}, eval(func)); if round(CC[method_cc]) = gear_cc then mt := 'gear' elif round(CC[method_cc]) = mgear_cc then mt := 'mgear' elif round(CC[method_cc]) = dverk78_cc then mt := 'dverk78' else mt := 'rkf45' fi; ff := dsolve({seq(Diff(v[i],v[1]) = ff[i-1],i= 2..m+1), seq(v[i](init[1]) = init[i], i = 2..m+1)}, {seq(v[i](v[1]),i= 2..m+1)}, type = numeric, output = listprocedure, method = mt); ff := subs(ff,[v[1], seq(v[j](v[1]),j= 2..m+1)]); t := init[1] - h; nb := 0; nf := 0; while t >= CC[ti_cc] do yy := traperror(ff(t)); if lasterror = yy then t := CC[ti_cc] - 1 else nb := nb - 1; B := [yy, op(B)] fi; t := t - h od; t := init[1] + h; while t <= CC[tf_cc] do yy := traperror(ff(t)); if lasterror = yy then t := CC[tf_cc] + 1 else nf := nf + 1; B := [op(B), yy] fi; t := t + h od; array(nb .. nf, B); end: rungekutta := proc(F) local h, m, B, CC, optargs, initvars, init, tt, vars, func, f, i, precision, t; options `Copyright 1993 by Daniel Schwalbe`; precision := max(Digits, trunc(evalhf(Digits))); CC := array(1..48); optargs := [seq(args[i], i = 2..nargs)]; if type(args[3],numeric) then t := op(1,[op(1,eval(F))]); optargs := ['inits' = args[2], 'stepsize' = args[3], 'numsteps' = args[4], t = args[2][1] .. args[2][1] + args[3]*args[4]]; if nargs > 4 then optargs := [op(optargs), 'iterations' = args[5]] fi; fi; optargs := map( proc(xx) if type(xx,equation) then xx else 'dummy' = xx fi end, optargs); initvars := `ODE/setup`(F, CC, optargs); vars := initvars[2]; init := initvars[1][1]; func := initvars[3]; `ODE/graphics/setcc`(F, CC, optargs, vars); f := `plots/BESIRK/makeproc`(func); B := `ODE/rkhfbandf`(f, CC, init); eval(B); end: makelist := proc(A, c10, c20) local c1, c2, ri, rf, i; options `Copyright 1991 by Daniel Schwalbe`; if nargs = 1 then ri := op(1, op(2, eval(A))); rf := op(2, op(2, eval(A))); c1 := 1; c2 := 2; else c1 := c10; c2 := c20; ri := op(1, op(2, eval(A))); rf := op(2, op(2, eval(A))) fi; [seq([A[i][c1], A[i][c2]], i = ri..rf)]; end: ccomputelineshf := proc(f, Ax, Ay, CC, m, n) local dx, dy, slp, xinc, yinc, pt, p, q, r, s, i, j, asp, boxx, boxy, blocks; p := CC[p_cc]; q := CC[q_cc]; r := CC[r_cc]; s := CC[s_cc]; dx := (q - p)/m; dy := (s - r)/(n); blocks := max(round(m/sqrt(2)), round(sqrt(2)*n)); asp := CC[asp_cc]; if asp < 1 then boxx := (q - p)/blocks*asp; boxy := (s - r)/blocks; else boxx := (q - p)/blocks; boxy := (s - r)/blocks/asp; fi; pt := array(1..2); slp := array(1..1); pt[1] := p + dx/2; for i from 0 to m-1 do pt[2] := r + (modp(i, 2) + 1)*dy/2; for j from 0 to n - 1 do f(pt, slp); xinc := `ODE/gfish`(boxx, boxy, 1, slp[1])/3; yinc := xinc*slp[1]; Ax[i, j, 1] := pt[1] - xinc; Ay[i, j, 1] := pt[2] - yinc; Ax[i, j, 2] := pt[1] + xinc; Ay[i, j, 2] := pt[2] + yinc; pt[2] := pt[2] + dy od; pt[1] := pt[1] + dx od; end: ccomputelines := proc(f, Ax, Ay, p, q, r, s, m, n) local dx, dy, ldx, slp, ldy, xinc, yinc, a, b, i, j, dydx; dx := evalf((q - p)/m); dy := evalf((s - r)/(n+1)); dydx := evalf(dy/dx); ldx := evalf(dx/3); ldy := evalf(dy/3); a := p + dx/2; for i from 0 to m - 1 do b := r + (modp(i,2) + 1)*dy/2; for j from 0 to n - 1 do slp := traperror(evalf(f(a, b))); if lasterror = slp then xinc := 0; yinc := 0; elif not type(slp, 'numeric') then xinc := 0; yinc := 0; elif slp > 0 then if slp > dydx then xinc := ldy/slp; yinc := ldy else xinc := ldx; yinc := ldx*slp fi; else if slp < -dydx then xinc := -ldy/slp; yinc := -ldy else xinc := ldx; yinc := ldx*slp fi; fi; Ax[i, j, 1] := a - xinc; Ay[i, j, 1] := b - yinc; Ax[i, j, 2] := a + xinc; Ay[i, j, 2] := b + yinc; b := b + dy od; a := a + dx od; end: `ODE/lineshell` := proc(f, F, CC, optargs) local mm, nn, flines, Ax, Ay, ff, p, q, r, s, fff, fl, j, i, fcol, fc, fthick; mm := CC[hboxes_cc]; nn := CC[vboxes_cc]; mm := 2*round(round(sqrt(2)*mm)/2); nn := round(nn/sqrt(2)); CC[hboxes_cc] := mm; CC[vboxes_cc] := nn; p := CC[p_cc]; q := CC[q_cc]; r := CC[r_cc]; s := CC[s_cc]; if mm > 0 and nn > 0 and CC[field_cc] = 1 then fthick := round(CC[fthick_cc]); fcol := subs(select(xx -> (op(1,xx) = 'fieldcolor'),optargs), 'fieldcolor'); if fcol = 'fieldcolor' then fc := unapply(COLOUR(RGB, 0,0,0), xdum, ydum) elif (fcol = 'grays' or fcol = 'graydirection') then fc := proc(y) if y > 0 then COLOUR(RGB, 0.7,0.7,0.7) else COLOUR(RGB, 0.3,0.3,0.3) fi end; elif (fcol = 'colors' or fcol = 'colordirection') then fc := proc(xdum) if xdum > 0 then COLOUR(HUE, 0.7) else COLOUR(HUE, 0.3) fi end; elif not type(fcol, procedure) then fc := unapply(fcol, xdum) else fc := fcol fi; Ax := array(0..mm-1, 0..nn-1, 1..2); Ay := array(0..mm-1, 0..nn-1, 1..2); fl := traperror(evalhf(ccomputelineshf( f, var(Ax), var(Ay), var(CC), mm, nn))): if lasterror = fl then if type(eval(F), {procedure, numeric}) then fff := subs('_body' = F(op(1, eval(F))), '_arg' = op(1, eval(F)), proc(_arg) _body end): else fff := subs('_body' = subs(op(1,lhs(F)) = _AAy, op(2,lhs(F)) = _AAt, rhs(F)), '_arg' = op([_AAt,_AAy]), proc(_arg) _body end ): fi; ccomputelines(fff, Ax, Ay, p, q, r, s, mm, nn); fi; flines := seq(seq(CURVES( [[Ax[i, j, 1], Ay[i, j, 1]], [Ax[i, j, 2], Ay[i, j, 2]]] , fc(Ay[i, j, 2] - Ay[i, j, 1]), THICKNESS(fthick) ), i = 0..mm-1), j = 0..nn-1): else flines := NULL fi: flines end: directionfield := proc(F) local init, i, j, f, pra, CC, p, q, r, s, flows, precision, flines, initlist, backg, fcol, func, nclines, jac, vars, vv, tt, tl, xl, optargs, ni, B, nnn, initvars: options `Copyright 1991 by Daniel Schwalbe`; precision := max(Digits, trunc(evalhf(Digits))); optargs := map( proc(xx) if type(xx,equation) then xx else 'dummy' = xx fi end, [args[2..nargs]]); CC := array(1..48): CC[m_cc] := 1; initvars := `ODE/setup`(F, CC, optargs); vars := initvars[2]; initlist := initvars[1]; func := initvars[3]; optargs := ['view' = [vars[1], vars[2]], op(optargs)]; pra := subs(optargs, [vars[1], vars[2]]); if not type(pra, [range, range]) then tt := select(xx -> type(rhs(xx), range), optargs); if nops(tt) < 2 then ERROR(`must supply ranges for plotting`) else optargs := [vars[1] = rhs(tt[1]), vars[2] = rhs(tt[2]), op(optargs)] fi fi; CC[field_cc] := 1; vv := `ODE/graphics/setcc`(F, CC, optargs, vars)[1]; f := `plots/BESIRK/makeproc`(func); if round(CC[method_cc]) = 1 then jac := `plots/ODE/makejac`(func) else jac := 1 fi; flines := `ODE/lineshell`(f, F, CC, optargs); flows := NULL; i := 0; ni := nops(initlist); for init in initlist do i := i+1; if round(CC[intmethod_cc]) = 1 then B := besirkbandf(f, jac, CC, init); elif round(CC[intmethod_cc]) = 2 then B := `ODE/dsolve`(func, CC, init); else B := `ODE/rkhfbandf`(f, CC, init); fi; flows := flows, `ODE/makeflow`(B, CC, i, ni, vv, vars, optargs); od; p := CC[p_cc]; q := CC[q_cc]; r := CC[r_cc]; s := CC[s_cc]; if round(CC[nclines_cc]) = 1 then nnn := plots[implicitplot]( convert(subs(_AA[1] = xx, _AA[2] = yy, eval(func)),set), xx = p..q, yy = r .. s): else nnn := NULL: fi; backg := subs(select(xx -> (op(1,xx) = 'background'),optargs), background); if backg = 'background' then backg := NULL else backg := POLYGONS([[p,s],[p,r],[q,r],[q,s]], backg) fi; tl := convert(vv[1][1], string); xl := convert(vv[1][2], string); PLOT(VIEW(p..q, r..s), AXESLABELS(tl, xl), flows, flines, CURVES(op(op(1,nnn)), COLOUR(RGB,0.75,0.75,0.75), THICKNESS(2)), CURVES(op(op(2,nnn)), COLOUR(RGB,0.75,0.75,0.75), THICKNESS(2)), backg, AXESSTYLE(BOX), STYLE(PATCHNOGRID)): end: `ODE/gfish` := proc(a, b, c, d) if (c = 0 and d = 0) then 10^7 else a*b / sqrt((b*c)^2 + (a*d)^2) fi end: `ODE/gpfish` := proc(a, b, c, d) if (c = 0 and d = 0) then 10^7 else sqrt( a*b / sqrt((b*c)^2 + (a*d)^2))/2 fi end: `ODE/computevects` := proc(Ax, Ay, f, CC) local i, j, k, minlen, ap, bp, boxx, boxy, t, x, y, flog, m, n, blocks, p, q, r, s, fx, fy, asp, fhead, gg, ggp, fishs, fishl, tx, ty, ux, uy, pux, puy, xin, yout; options `Copyright 1991 by Daniel Schwalbe`; xin := array(1..3); yout := array(1..2); p := CC[p_cc]; q := CC[q_cc]; r := CC[r_cc]; s := CC[s_cc]; m := CC[hboxes_cc]; n := CC[vboxes_cc]; fhead := CC[fhead_cc]; fishl := CC[fishl_cc]; fishs := CC[fishs_cc]; blocks := max(m, n); minlen := 10^7; asp := CC[asp_cc]; if asp < 1 then boxx := (q - p)/blocks*asp; boxy := (s - r)/blocks; else boxx := (q - p)/blocks; boxy := (s - r)/blocks/asp; fi; for j from 0 to m do for i from 0 to n do xin[2] := Ax[j, 5*i]; xin[3] := Ay[j, 5*i]; f(xin, yout); Ax[j, 5*i + 1] := yout[1]; Ay[j, 5*i + 1] := yout[2]; gg := `ODE/gfish`(boxx, boxy, yout[1], yout[2]); Ax[j, 5*i + 2] := gg; if gg = 0 then gg := 10^7 fi; minlen := min(gg, minlen); od; od; fishl := evalf(log[10](max(min(fishl,10.),1.))); fishs := 4*fishs*minlen^(1-fishl); for j from 0 to m do for i from 0 to n do k := 5*i; tx := Ax[j, k + 1]; ty := Ay[j, k + 1]; gg := Ax[j, k + 2];; flog := evalf(fishs*exp(log(gg)*fishl)); ap := Ax[j, k]; bp := Ay[j, k]; ux := tx*flog/4; uy := ty*flog/4; pux := -fhead*ty*boxx/boxy*flog/4; puy := fhead*tx*boxy/boxx*flog/4; Ax[j, k + 1] := ap + ux; Ay[j, k + 1] := bp + uy; Ax[j, k + 2] := ap + ux*0.75 - 0.25*pux; Ay[j, k + 2] := bp + uy*0.75 - 0.25*puy; Ax[j, k + 3] := ap + ux*0.75 + 0.25*pux; Ay[j, k + 3] := bp + uy*0.75 + 0.25*puy; Ax[j, k + 4] := 1 - minlen/gg; od;od end: `ODE/computefish` := proc(Ax, Ay, f, df, CC) local i, j, k, minlen, ap, bp, boxx, boxy, t, x, y, flog, m, n, gg2, blocks, dx, dy, p, q, r, s, f2x, f2y, f3x, f3y, txp, typ, asp, fhead, gg, ggp, fishs, fishl, tx, ty, ux, uy, pux, puy, xin, yout; options `Copyright 1991 by Daniel Schwalbe`; xin := array(1..3); yout := array(1..2); p := CC[p_cc]; q := CC[q_cc]; r := CC[r_cc]; s := CC[s_cc]; m := CC[hboxes_cc]; n := CC[vboxes_cc]; fhead := CC[fhead_cc]/2; fishl := CC[fishl_cc]; fishs := CC[fishs_cc]; dx := (q - p)/m: dy := (s - r)/n: blocks := max(m, n); minlen := 10^7; asp := CC[asp_cc]; if asp < 1 then boxx := (q - p)/blocks*asp; boxy := (s - r)/blocks; else boxx := (q - p)/blocks; boxy := (s - r)/blocks/asp; fi; for j from 0 to m do for i from 0 to n do xin[2] := Ax[j, 9*i]; xin[3] := Ay[j, 9*i]; f(xin, yout); Ax[j, 9*i + 1] := yout[1]; Ay[j, 9*i + 1] := yout[2]; gg := `ODE/gfish`(boxx, boxy, yout[1], yout[2]); Ax[j, 9*i + 2] := gg; if gg = 0 then gg := 10^7 fi; minlen := min(gg, minlen); od; od; fishl := evalf(log[10](max(min(fishl,10.),1.))); fishs := 2*fishs*minlen^(1-fishl) ; for j from 0 to m do for i from 0 to n do k := 9*i; xin[2] := Ax[j, k]; xin[3] := Ay[j, k]; df(xin, yout); txp := yout[1]; typ := yout[2]; tx := Ax[j, 9*i + 1]; ty := Ay[j, 9*i + 1]; gg := Ax[j, 9*i + 2]; ux := fhead*tx*gg/8; uy := fhead*ty*gg/8; pux := -fhead*ty*boxx/boxy*gg/8; puy := fhead*tx*boxy/boxx*gg/8; ggp := min(gg, `ODE/gpfish`(boxx, boxy, txp, typ)); flog := evalf(fishs*exp(log(ggp)*fishl)) ; ap := Ax[j, k]; bp := Ay[j, k]; f2x := tx*flog/2 - (flog/2)^2/2*txp; f2y := ty*flog/2 - (flog/2)^2/2*typ; f3x := tx*flog - (flog)^2/2*txp; f3y := ty*flog - (flog)^2/2*typ; Ax[j, k] := ap + ux; Ay[j, k] := bp + uy; Ax[j, k + 1] := ap + .707*ux + .707*pux; Ay[j, k + 1] := bp + .707*uy + .707*puy; Ax[j, k + 2] := ap + pux; Ay[j, k + 2] := bp + puy; Ax[j, k + 3] := ap - f2x + pux/2; Ay[j, k + 3] := bp - f2y + puy/2; Ax[j, k + 4] := ap - f3x; Ay[j, k + 4] := bp - f3y; Ax[j, k + 5] := ap - f2x - pux/2; Ay[j, k + 5] := bp - f2y - puy/2; Ax[j, k + 6] := ap - pux; Ay[j, k + 6] := bp - puy; Ax[j, k + 7] := ap + .707*ux - .707*pux; Ay[j, k + 7] := bp + .707*uy - .707*puy; Ax[j, k + 8] := 1 - minlen/gg; od;od end: phasevectsShell := proc(f, df, CC, optargs) local Ax, Ay, pv, mm, nn, phasevects,aa, j, i, p, q, r, s, fc, fcol, fthick; mm := round(CC[hboxes_cc]); nn := round(CC[vboxes_cc]); p := CC[p_cc]; q := CC[q_cc]; r := CC[r_cc]; s := CC[s_cc]; fthick := round(CC[fthick_cc]); fcol := subs(select(xx -> (op(1,xx) = 'fieldcolor'),optargs), 'fieldcolor'); if fcol = 'fieldcolor' then fc := unapply(COLOUR(RGB, 0,0,0), xdum); elif not type(fcol,procedure) then fc := unapply(fcol,xdum) else fc := fcol fi; if CC[field_cc] = 2 then mm := 2*round(round(sqrt(2)*mm)/2); nn := round(nn/sqrt(2)); CC[hboxes_cc] := mm; CC[vboxes_cc] := nn; Ax := array(0..mm, 0..5*nn + 4); Ay := array(0..mm, 0..5*nn + 4); for i from 0 to nn do for j from 0 to round(mm/2) do Ax[2*j, 5*i] := evalf(p + (0.5 + i)*(q - p)/(nn+1.5)); Ay[2*j, 5*i] := evalf(r + (1 + 2*j)*(s - r)/(mm+1.5)); od od; for i from 0 to nn do for j from 0 to round(mm/2)-1 do Ax[2*j+1, 5*i] := evalf(p + (1 + i)*(q - p)/(nn+1.5)); Ay[2*j+1, 5*i] := evalf(r + (2 + 2*j)*(s - r)/(mm+1.5)); od od; pv := traperror(evalhf(`ODE/computevects`(var(Ax), var(Ay), f, var(CC)))); if lasterror = pv then `ODE/computevects`(Ax,Ay,f,CC); fi; phasevects := seq(seq( CURVES( [[Ax[j, 5*i], Ay[j, 5*i]], [Ax[j, 5*i + 1], Ay[j, 5*i + 1]], [Ax[j, 5*i + 2], Ay[j, 5*i + 2]], [Ax[j, 5*i + 1], Ay[j, 5*i + 1]], [Ax[j, 5*i + 3], Ay[j, 5*i + 3]]], fc(Ax[j, 5*i + 4]/2), THICKNESS(fthick)), j = 0..mm), i = 0..nn) else Ax := array(0..mm, 0..9*nn + 8); Ay := array(0..mm, 0..9*nn + 8); for j from 0 to mm do for i from 0 to nn do Ax[j, 9*i] := evalf(p + rand()/10^12*(q - p)); Ay[j, 9*i] := evalf(r + rand()/10^12*(s - r)) od od; pv := traperror(evalhf(`ODE/computefish`(var(Ax), var(Ay), f, df, var(CC)))); if lasterror = pv then `ODE/computefish`(Ax,Ay,f,df,CC); fi; phasevects := seq(seq(op(subs(aa=fc(Ax[j, 9*i + 8]/2), [ POLYGONS( [[Ax[j, 9*i], Ay[j, 9*i]], [Ax[j, 9*i + 1], Ay[j, 9*i + 1]], [Ax[j, 9*i + 2], Ay[j, 9*i + 2]], [Ax[j, 9*i + 3], Ay[j, 9*i + 3]], [Ax[j, 9*i + 5], Ay[j, 9*i + 5]], [Ax[j, 9*i + 6], Ay[j, 9*i + 6]], [Ax[j, 9*i + 7], Ay[j, 9*i + 7]]], aa), POLYGONS( [[Ax[j, 9*i + 3], Ay[j, 9*i + 3]], [Ax[j, 9*i + 4], Ay[j, 9*i + 4]], [Ax[j, 9*i + 5], Ay[j, 9*i + 5]]], aa)])), j = 0..mm), i = 0..nn) fi; phasevects end: `ODE/makeflow` := proc(B, CC, i, ni, viewt, vars, optargs) local ii, j, k, l, Arkx, Arky, nb, nf, BB, BBB, fc, xdum, n, npts, othick, v1, v2, v3, m, vvv, vv, vc, fhead, ti, tf, viewnum, fcol, p, q, hm, ccol; m := CC[m_cc]; ti := CC[ti_cc]; tf := CC[tf_cc]; fhead := CC[fhead_cc]; p := CC[p_cc]; q := CC[q_cc]; othick := round(CC[orbthick_cc]); fcol := subs(select(xx -> (op(1,xx) = 'orbitcolor'),optargs), 'orbitcolor'); if fcol = 'orbitcolor' then fcol := subs(select(xx -> (op(1, xx) = 'flowcolor'),optargs), 'flowcolor'); fi; if (fcol = 'flowcolor' or fcol = 'orbitcolor') then fc := xdum -> COLOUR(RGB, 0,0,0); else if fcol = 'RAINBOW' or fcol = 'HUE' then fc := unapply(COLOUR(HUE, (1+xdum)*0.45), xdum) elif fcol = 'GRAY' then fc := unapply(COLOUR(RGB, (1-xdum)*0.75+0.01, (1-xdum)*0.75+0.01, (1-xdum)*0.75+0.01), xdum); elif not type(fcol,procedure) then fc := unapply(fcol,xdum) else fc := fcol fi; fi; nf := op(2,op(2,eval(B))); nb := op(1,op(2,eval(B))); viewnum := nops(viewt); BB := NULL; vc := 0; for vv in viewt do vc := vc + 1; vvv := []; for ii to nops(vv) do vvv := [op(vvv), op(select( unapply(vars[xx] = vv[ii],xx) , [seq(j,j=1..m+1)]))]; od; v1 := vvv[1]; v2 := vvv[2]; if nops(vvv) = 3 then v3 := vvv[3] fi; if nops(vvv) = 2 then if round(CC[flowf_cc]) = 2 then hm := fhead*(q-p)/40/CC[hmax_cc]/8; n := trunc((nf-nb)/CC[segments_cc]); if n > 0 then npts := round(CC[segments_cc]); if v1 = 2 then BB := BB, seq(op([CURVES( [[B[(j+1)*npts+nb][2] - (B[(j+1)*npts+nb][2] - B[(j+1)*npts+nb-1][2])*hm, (B[(j+1)*npts+nb][2] - B[(j+1)*npts+nb-1][2])*hm], [B[(j+1)*npts+nb][2], 0], [B[(j+1)*npts+nb][2]- (B[(j+1)*npts+nb][2] - B[(j+1)*npts+nb-1][2])*hm, -(B[(j+1)*npts+nb][2] - B[(j+1)*npts+nb-1][2])*hm]], fc((abs(B[j*npts + 1 + nb][2]-B[(j+1)*npts + nb][2]))/(tf-ti)), THICKNESS(othick)), CURVES([seq( [B[j*npts + ii+nb][2], 0], ii = 1 .. npts)], fc((abs(B[j*npts + 1 + nb][2]-B[(j+1)*npts + nb][2]))/(tf-ti)), THICKNESS(othick)) ]), j = 0..n-1) else BB := BB, seq(op([CURVES( [[(B[(j+1)*npts+nb][2] - B[(j+1)*npts+nb-1][2])*hm, B[(j+1)*npts+nb][2] - (B[(j+1)*npts+nb][2] - B[(j+1)*npts+nb-1][2])*hm], [0, B[(j+1)*npts+nb][2]], [-(B[(j+1)*npts+nb][2] - B[(j+1)*npts+nb-1][2])*hm, B[(j+1)*npts+nb][2] - (B[(j+1)*npts+nb][2] - B[(j+1)*npts+nb-1][2])*hm]], fc(ccol), THICKNESS(othick)), CURVES([seq( [0, B[j*npts + ii+nb][2]], ii = 1 .. npts) ],fc(ccol), THICKNESS(othick))]), j = 0..n-1) fi fi elif CC[flowf_cc] = 0 then if fcol = 'HUE' or fcol = 'GRAY' then BB := BB, seq(CURVES([[B[j][v1],B[j][v2]],[B[j+1][v1],B[j+1][v2]]], fc((B[j][1]-ti)/(tf-ti)), THICKNESS(othick)), j = nb .. nf-1); elif fcol = 'RAINBOW' then BB := BB, CURVES( [seq([B[j][v1],B[j][v2]], j = nb .. nf)], fc((viewnum*(i-1) + vc)/(viewnum*ni)), THICKNESS(othick)); else BB := BB, CURVES( [seq([B[j][v1],B[j][v2]], j = nb .. nf)], fc((viewnum*(i-1) + vc)/(viewnum*ni)), THICKNESS(othick)); fi; else n := trunc((nf-nb)/CC[segments_cc]); if n > 0 then npts := round(CC[segments_cc]); Arkx := array(0..n-1, 1..npts+1,1..2); Arky := array(0..n-1, 1..npts+1,1..2); BBB := array(nb .. nf, 2..3); for k from nb to nf do BBB[k,2] := B[k][v1] od; for k from nb to nf do BBB[k,3] := B[k][v2] od; evalhf(parflowlist(var(BBB), var(Arkx), var(Arky), var(CC), nb, nf)); BB := BB, seq(op([POLYGONS([ [Arkx[j,npts+1,1] - Arkx[j,npts+1,2], Arky[j,npts+1,1]], [Arkx[j,npts+1,1] - 0.707*Arkx[j,npts+1,2], Arky[j,npts+1,1] + 0.707*Arky[j,npts+1,2]], [Arkx[j,npts+1,1], Arky[j,npts+1,1] + Arky[j,npts+1,2]], [Arkx[j,npts+1,1] + 0.707*Arkx[j,npts+1,2], Arky[j,npts+1,1] + 0.707*Arky[j,npts+1,2]], [Arkx[j,npts+1,1] + Arkx[j,npts+1,2], Arky[j,npts+1,1]], [Arkx[j,npts+1,1] + 0.707*Arkx[j,npts+1,2], Arky[j,npts+1,1] - 0.707*Arky[j,npts+1,2]], [Arkx[j,npts+1,1], Arky[j,npts+1,1] - Arky[j,npts+1,2]], [Arkx[j,npts+1,1] - 0.707*Arkx[j,npts+1,2], Arky[j,npts+1,1] - 0.707*Arky[j,npts+1,2]] ],fc(1)), seq(POLYGONS([ [Arkx[j,ii,1], Arky[j,ii,1]], [Arkx[j,ii,2], Arky[j,ii,2]], [Arkx[j,ii+1,2], Arky[j,ii+1,2]], [Arkx[j,ii+1,1], Arky[j,ii+1,1]] ],fc((ii)/npts)), ii = 1 .. npts-1)]), j = 0..n-1) fi fi else if fcol <> 'HUE' and fcol <> 'GRAY' then BB := BB, CURVELIST( [seq([B[j][v1],B[j][v2],B[j][v3]], j = nb .. nf)], fc((viewnum*(i-1) + vc)/(viewnum*ni)), THICKNESS(othick)); else BB := BB, MESH( [seq([[B[j][v1],B[j][v2],B[j][v3]], [B[j][v1],B[j][v2],B[j][v3]]], j = nb .. nf)]); fi; fi; od; BB end: parflowlist := proc(Ain, Ax, Ay, CC, nb, nf) local j, n, npts, asp, fishx, fishy, fhead, boxx, boxy, nj, gg, dx, dy, p, q, r, s: options `Copyright 1995 by Daniel Schwalbe`; fhead := CC[fhead_cc]/8.; p := CC[p_cc]; q := CC[q_cc]; r := CC[r_cc]; s := CC[s_cc]; asp := CC[asp_cc]; if asp < 1 then boxx := (q - p)/20*asp; boxy := (s - r)/20; else boxx := (q - p)/20; boxy := (s - r)/20/asp; fi; n := trunc((nf-nb)/CC[segments_cc]); npts := round(CC[segments_cc]); fishx := array(1..npts,1..2); fishy := array(1..npts,1..2); for j from 0 to n-1 do for nj to npts do dx := Ain[j*npts + nj + nb,2] - Ain[j*npts + nj - 1 + nb,2]; dy := Ain[j*npts + nj + nb,3] - Ain[j*npts + nj - 1 + nb,3]; gg := `ODE/gfish`(boxx, boxy, dx, dy); fishx[nj,1] := Ain[j*npts + nj + nb,2]; fishy[nj,1] := Ain[j*npts + nj + nb,3]; fishx[nj,2] := -fhead*gg*boxx/boxy*dy*(nj-1)/(npts-1); fishy[nj,2] := fhead*gg*boxy/boxx*dx*(nj-1)/(npts-1) od: for nj to npts do Ax[j, nj, 1] := fishx[nj,1] + fishx[nj,2]; Ax[j, nj, 2] := fishx[nj,1] - fishx[nj,2]; Ay[j, nj, 1] := fishy[nj,1] + fishy[nj,2]; Ay[j, nj, 2] := fishy[nj,1] - fishy[nj,2]; od; Ax[j,npts+1,1] := fishx[npts,1]; Ay[j,npts+1,1] := fishy[npts,1]; Ax[j,npts+1,2] := -fhead*boxx; Ay[j,npts+1,2] := fhead*boxy; od: end: phaseline := proc(F) local B, CC, optargs, vv, ff, i, j, flows, backg, xl, signz, m, f, initlist, init, func, vars, initvars, ni, prange, vert, p, r, q, s, nf, nb, ddone, FF, axl, axt, zeroes, z1, z2, ti, tf, npts, xnf, xnb, flowtype, tt; optargs := map( proc(xx) if type(xx,equation) then xx else 'dummy' = xx fi end, [args[2..nargs]]); if type(F, procedure) then vars := [op(1,eval(F))]; elif type(F[1], procedure) then vars := [op(1,eval(F[1]))]; else if not type(F,list) then FF := [F] else FF := F fi; vars := subs(optargs, variables); if vars = 'variables' then vars := [op(2,lhs(FF[1])), op(1,lhs(FF[1]))]; fi; fi; vars := [vars[1], vars[2], 'ydum']; prange := subs(optargs, vars[2]); if not type(prange, range) then tt := select(xx -> type(rhs(xx), range) ,optargs); if nops(tt) > 0 then prange := rhs(tt[1]) else prange := -10 .. 10 fi fi; optargs := [aspectratio = 0.1, numsteps = 50, segments = 5, op(optargs), 'view' = [vars[2], vars[3]], 'bound' = true, 'vertical' = false, 'flowcolor' = black, inits = [0,1], 'tdummy' = -2.5..2.5]; CC := array(1..48): CC[p_cc] := op(1, prange); CC[q_cc] := op(2, prange); CC[s_cc] := (CC[q_cc] - CC[p_cc])/20; CC[r_cc] := -CC[s_cc]; p := CC[p_cc]; q := CC[q_cc]; r := CC[r_cc]; s := CC[s_cc]; CC[m_cc] := 1; initvars := `ODE/setup`(F, CC, optargs); initlist := initvars[1]; func := initvars[3]; func := array([func[1],0]); CC[m_cc] := 2; `ODE/graphics/setcc`(F, CC, optargs, vars); flowtype := subs(select(xx -> (op(1,xx) = 'flowfield'),optargs), flowfield); if flowtype = 'flowfield' then CC[flowf_cc] := 2 else CC[flowf_cc] := 1 fi; vert := subs(optargs,'vertical'); if vert then p := CC[r_cc]; q := CC[s_cc]; r := CC[p_cc]; s := CC[q_cc]; vv := [[vars[3], vars[2]]]; xl := convert(vv[1][2], string); axl := AXESLABELS(` `,xl); axt := AXESTICKS(0, DEFAULT); else vv := [[vars[2], vars[3]]]; xl := convert(vv[1][1], string); axl := AXESLABELS(xl, ` `); axt := AXESTICKS(DEFAULT, 0); fi; f := `plots/BESIRK/makeproc`(func); flows := NULL; i := 0; ni := 10; ddone := false; init := [0, CC[p_cc], 0]; CC[fhead_cc] := 8*CC[fhead_cc]; npts := round(CC[segments_cc]); CC[hmax_cc] := 0; ti := CC[ti_cc]; tf := CC[tf_cc]; CC[hmin_cc] := (CC[q_cc]-CC[p_cc])/80; zeroes := {fsolve(func[1], _AA[2])}; while not ddone do i := i+1; if evalf(subs(_AA[2] = init[2], func[1])) > 0 then CC[ti_cc] := 0; CC[tf_cc] := tf - ti else CC[ti_cc] := ti - tf; CC[tf_cc] := 0 fi; B[i] := `ODE/rkhfbandf`(f, CC, init); nf := op(2,op(2,eval(B[i]))); nb := op(1,op(2,eval(B[i]))); xnf := B[i][nf][2]; xnb := B[i][nb][2]; init := [0, max(xnf, xnb) + (CC[q_cc]-CC[p_cc])/40, 0]; z1 := fsolve(func[1], _AA[2], min(xnf, xnb) .. max(xnf,xnb)); z2 := fsolve(func[1], _AA[2], max(xnf, xnb) .. init[2]); zeroes := {op(zeroes), z1, z2}; if init[2] >= CC[q_cc] then ddone := true fi; od; for j to i do flows := flows, `ODE/makeflow`(B[j], CC, j, ni, vv, vars, optargs) od; zeroes := select(xx -> type(xx,numeric), zeroes); if vert then zeroes := sort(select( subs({a=r, b=s}, proc(xx) evalb(xx <= b and xx >= a) end), [r,op(zeroes),s])) else zeroes := sort(select( subs({a=p, b=q}, proc(xx) evalb(xx <= b and xx >= a) end), [p,op(zeroes),q])) fi; signz := seq(signum(evalf(subs(_AA[2] = (zeroes[i+1] + zeroes[i])/2, func[1]))),i = 1 .. nops(zeroes)-1); backg := subs(select(xx -> (op(1,xx) = 'background'),optargs), background); if backg = 'background' then backg := NULL; for i to nops(zeroes) - 1 do if vert then backg := backg, POLYGONS([[q,zeroes[i]], [p,zeroes[i]], [p,zeroes[i+1]], [q,zeroes[i+1]]], COLOR(RGB, (signz[i] + 4)/8,(signz[i] + 4)/8,(signz[i] + 4)/8)) else backg := backg, POLYGONS([[zeroes[i],s], [zeroes[i],r], [zeroes[i+1],r], [zeroes[i+1],s]], COLOR(RGB, (signz[i] + 4)/8,(signz[i] + 4)/8,(signz[i] + 4)/8)) fi od else backg := POLYGONS([[p,s],[p,r],[q,r],[q,s]], backg) fi; PLOT(VIEW(p..q, r..s), axl, flows, backg, AXESSTYLE(BOX), SCALING(CONSTRAINED), STYLE(PATCHNOGRID), axt) end: phaseplot := proc(F) local i, p, q, r, s, phasevects, initlist, nnn, init, B, f, optargs, ni, initvars, pra, CC, func, nclines, df, vars, flows, jac, backg, vv, tt, tl, xl: options `Copyright 1991 by Daniel Schwalbe`; optargs := map( proc(xx) if type(xx,equation) then xx else 'dummy' = xx fi end, [args[2..nargs]]); CC := array(1..48): CC[m_cc] := 2; CC[field_cc] := 2; optargs := [op(optargs), 'tdummy' = -2.5..2.5]; initvars := `ODE/setup`(F, CC, optargs); vars := initvars[2]; initlist := initvars[1]; func := initvars[3]; optargs := ['view' = [vars[2], vars[3]], op(optargs)]; pra := subs(optargs, [vars[2], vars[3]]); if not type(pra, [range, range]) then tt := select(xx -> type(rhs(xx), range), optargs); if nops(tt) < 2 then ERROR(`must supply ranges for plotting`) else optargs := [vars[2] = rhs(tt[1]), vars[3] = rhs(tt[2]), op(optargs)] fi fi; vv := `ODE/graphics/setcc`(F, CC, optargs, vars)[1]; f := `plots/BESIRK/makeproc`(func); if round(CC[method_cc]) = 1 then jac := `plots/ODE/makejac`(func) else jac := 1 fi; df := array([func[1]*diff(func[1], _AA[2]) + func[2]*diff(func[1], _AA[3]), func[1]*diff(func[2], _AA[2]) + func[2]*diff(func[2], _AA[3])]); df := `plots/BESIRK/makeproc`(df); if CC[hboxes_cc] > 0 and CC[vboxes_cc] > 0 and CC[field_cc] > 0 then phasevects := phasevectsShell(f, df, CC, optargs); else phasevects := NULL fi; flows := NULL; i := 0; ni := nops(initlist); for init in initlist do i := i+1; if round(CC[intmethod_cc]) = 1 then B := besirkbandf(f, jac, CC, init); elif round(CC[intmethod_cc]) = 2 then B := `ODE/dsolve`(func, CC, init); else B := `ODE/rkhfbandf`(f, CC, init); fi; flows := flows, `ODE/makeflow`(B, CC, i, ni, vv, vars, optargs); od; p := CC[p_cc]; q := CC[q_cc]; r := CC[r_cc]; s := CC[s_cc]; if round(CC[nclines_cc]) = 1 then nnn := plots[implicitplot]( convert(subs(_AA[2] = xx, _AA[3] = yy, eval(func)),set), xx = p..q, yy = r .. s): else nnn := NULL: fi; backg := subs(select(xx -> (op(1,xx) = 'background'),optargs), background); if backg = 'background' then backg := NULL else backg := POLYGONS([[p,s],[p,r],[q,r],[q,s]], backg) fi; tl := convert(vv[1][1], string); xl := convert(vv[1][2], string); PLOT(VIEW(p..q, r..s), AXESLABELS(tl, xl), flows, phasevects, CURVES(op(op(1,nnn)), COLOUR(RGB,0.75,0.75,0.75), THICKNESS(2)), CURVES(op(op(2,nnn)), COLOUR(RGB,0.75,0.75,0.75), THICKNESS(2)), backg, AXESSTYLE(BOX), STYLE(PATCHNOGRID) ): end: picard := proc(f, init, n) local parray, ptemp, t0, y0: options `Copyright 1991 by Daniel Schwalbe`; t0 := init[1]; y0 := init[2]; ptemp := y0 + int(f('s', y0), 's' = t0..'t'); parray := ptemp; from 2 to n do ptemp := y0 + int(f('s', subs('t' = 's', ptemp)), 's' = t0..'t'); parray := parray, sort(ptemp) od; parray; end: `ODE/getnumdiff` := proc(F, optargs, m) local numdiff, FF, i, vars; if type(F, procedure) then numdiff := m; elif type(F[1], procedure) then numdiff := m; else if not type(F,list) then FF := [F] else FF := F fi; vars := subs(optargs,variables); if vars <> 'variables' then numdiff := 0; for i to m do if op(0,(FF[i])) = diff or op(0,lhs(FF[i])) = Diff then numdiff := numdiff + 1; fi; od; else numdiff := m; fi: fi; numdiff; end: `plots/BESIRK/makeproc` := proc(f) local i, j, m; m := op(2, op(2, eval(f))); readlib(procmake)( `&proc`([_AA, _eq], [], [], `&statseq`(seq( `&:=`(`&args`[2][i], f[i]), i = 1..m)) )); end: `plots/ODE/makejac` := proc(f) local i, j, k, jac, m; m := op(2, op(2, eval(f))); jac := linalg[jacobian](f, [seq(_AA[j], j = 2..m + 1)]); readlib(procmake)( `&proc`([_AA, _eq], [], [], `&statseq`(seq(seq( `&:=`(`&args`[2][i, j], jac[i, j]), i = 1..m), j = 1..m)) )) end: ludecomphf := proc (A,indx,n) local i,j,k,tiny,_A,S,Amax,VV,d,dum,imax; tiny := 1e-20; d := 1; VV := array(1..n); for i from 1 to n do Amax := 0; for j from 1 to n do if abs(A[i,j]) > Amax then Amax := abs(A[i,j]) fi; od; if Amax = 0 then ERROR (`Simgular Matrix`) fi; VV[i] := 1 / Amax; od; for j from 1 to n do for i from 1 to j-1 do S := A[i,j]; for k from 1 to i-1 do S := S - A[i,k] * A[k,j] od; A[i,j] := S; od; Amax := 0; for i from j to n do S := A[i,j]; for k from 1 to j-1 do S := S - A[i,k] * A[k,j]; od; A[i,j] := S; dum := VV[i] * abs(S); if dum >= Amax then imax := i; Amax := dum; fi; od; if j <> imax then for k from 1 to n do dum := A[imax,k]; A[imax,k] := A[j,k]; A[j,k] := dum; od; d = - d; VV[imax] := VV[j]; fi; indx[j] := imax; if A[j,j] = 0 then A[j,j] := tiny; fi; if j <> n then dum := 1 / A[j,j]; for i from j+1 to n do A[i,j] := A[i,j] * dum; od; fi; od; end: lubacksubhf := proc(A, b, indx, n) local i,j,k,ii,ll,S; ii := 0; for i from 1 to n do ll := indx[i]; S := b[ll]; b[ll] := b[i]; if ii <> 0 then for j from ii to i-1 do S := S - A[i,j] * b[j]; od; elif S <> 0 then ii := 1; fi; b[i] := S; od; for i from n to 1 by -1 do S := b[i]; if i < n then for j from i+1 to n do S := S - A[i,j] * b[j]; od; fi; b[i] := S / A[i,i]; od; end: `BESIRK/polpol` := proc(CC, DidX, DidY, X, YY, dY) local i, j, m, suse, s, m1, k1, delta, q, f1, f2, Dw; m := CC[m_cc]; suse := CC[sUse_cc]; s := CC[ss_cc]; Dw := array(1..m); DidX[s] := X; for j to m do dY[j] := YY[CC[laststored_cc], j] od; if s = 1 then for j to m do DidY[j, 1] := YY[CC[laststored_cc], j] od else if s < suse then m1 := s else m1 := suse fi; for j to m do Dw[j] := YY[CC[laststored_cc], j] od; for k1 to m1 - 1 do delta := 1/(DidX[s - k1] - X); f1 := X*delta; f2 := DidX[s - k1]*delta; for j to m do q := DidY[j, k1]; DidY[j, k1] := dY[j]; delta := Dw[j] - q; if abs(delta) > 10^20 then CC[Idid_cc] := polerr; else ; dY[j] := f1*delta; Dw[j] := f2*delta; YY[CC[laststored_cc], j] := YY[CC[laststored_cc], j] + dY[j];fi od; if CC[Idid_cc] <> polerr then for j to m do DidY[j, s] := dY[j] od; fi; od; fi; end: StandT := proc(CC, YY, dY) local AbsTol, RelTol, m, t, et, finished, i, o, l; m := CC[m_cc]; AbsTol := CC[AbsTol_cc]; RelTol := CC[RelTol_cc]; finished := true; et := 10^(-30); for i to m do o := abs(dY[i]); l := AbsTol + RelTol*abs(YY[CC[laststored_cc], i]) + 10^(-30); et := max(et, o/l); if o > l then finished := false fi; od; CC[Ecur_cc] := et; finished; end: `BESIRK/rk4` := proc(f, YY, CC) local j, i, h, n, npts, Ytemp, k1, k2, k3, k4, m, y: options `Copyright 1993 by Daniel Schwalbe`; h := CC[h_cc]; CC[Idid_cc] := 0; npts := `besirk/sSeq`[CC[sUsed_cc]]; m := CC[m_cc]; Ytemp := array(1..m + 1); y := array(1..m + 1); for i to m + 1 do Ytemp[i] := YY[CC[laststored_cc], i - 1] od; k1 := array(1..m);k2 := array(1..m); k3 := array(1..m);k4 := array(1..m); to npts do f(Ytemp, k1): y[1] := Ytemp[1] + h/2; for i from 2 to m + 1 do y[i] := Ytemp[i] + h*k1[i - 1]/2 od; f(y, k2): for i from 2 to m + 1 do y[i] := Ytemp[i] + h*k2[i - 1]/2 od; f(y, k3): y[1] := Ytemp[1] + h; for i from 2 to m + 1 do y[i] := Ytemp[i] + h*k3[i - 1] od; f(y, k4): Ytemp[1] := y[1]; for i from 2 to m + 1 do Ytemp[i] := Ytemp[i] + h*(k1[i - 1] + 2*k2[i - 1] + 2*k3[i - 1] + k4[i - 1])/6 od: od: CC[laststored_cc] := CC[laststored_cc] + 1; for i from 0 to m do YY[CC[laststored_cc], i] := Ytemp[i + 1]; od end: `BESIRK/rk2` := proc(f, YY, CC) local j, i, h, n, npts, Ytemp, k1, k2, m, y: options `Copyright 1993 by Daniel Schwalbe`; h := CC[h_cc]; CC[Idid_cc] := 0; npts := `besirk/sSeq`[CC[sUsed_cc]]; m := CC[m_cc]; Ytemp := array(1..m + 1); y := array(1..m + 1); for i to m + 1 do Ytemp[i] := YY[CC[laststored_cc], i - 1] od; k1 := array(1..m);k2 := array(1..m); to npts do f(Ytemp, k1): y[1] := Ytemp[1] + h; for i from 2 to m + 1 do y[i] := Ytemp[i] + h*k1[i - 1] od; f(y, k2): Ytemp[1] := y[1]; for i from 2 to m + 1 do Ytemp[i] := Ytemp[i] + h*(k1[i - 1] + k2[i - 1])/2 od: od: CC[laststored_cc] := CC[laststored_cc] + 1; for i from 0 to m do YY[CC[laststored_cc], i] := Ytemp[i + 1] ; od end: `BESIRK/msirk3`:=proc(f,jac,YY,CC,sSeq) local a,b2,c2,b31,b32,R1,R2,R3,Ytemp, Y,BhaJ,m,n,h,i,j,k,Jac,F,k1,k2,k3,indx; m:=CC[m_cc]; h:=CC[h_cc]; CC[Idid_cc]:=0; n := `besirk/sSeq`[CC[sUsed_cc]]; a := CC[a_cc]; b2 := 3/4; c2 := 27/32; b31 := -(1/6) * (-2*a+1+8*a^2)/a; b32 := (2/9) *(-6*a+6*a^2+1)/a; R1 := 11/27 - b31; R2 := 16/27 - b32; R3 := 1; Y:=array(1..m+1); Ytemp:=array(1..m+1); Jac:=array(1..m,1..m); F:=array(1..m); BhaJ:=array(1..m,1..m); k1 := array(1..m); k2 := array(1..m); k3 := array(1..m); indx := array(1..m); for i to m+1 do Ytemp[i]:=YY[CC[laststored_cc],i-1] od; for k to n do for i to m+1 do Y[i]:=Ytemp[i] od; jac(Y,Jac); for i to m do for j to m do if (i=j and i<=CC[numdiff_cc]) then BhaJ[i,j]:=1-a*h*Jac[i,j] else BhaJ[i,j]:=-a*h*Jac[i,j] fi; od od; ludecomphf(BhaJ,indx,m); f(Y,F); for i to m do k1[i] := h * F[i]; od; lubacksubhf(BhaJ,k1,indx,m); Y[1]:=Y[1]+c2*h; for i to m do Y[i+1]:=Y[i+1]+b2*k1[i] od; f(Y,F); for i to m do k2[i] := h * F[i]; od; lubacksubhf(BhaJ,k2,indx,m); for i to CC[numdiff_cc] do k3[i] := b31 * k1[i] + b32 * k2[i]; od; for i from CC[numdiff_cc]+1 to m do k3[i]:=0 od; lubacksubhf(BhaJ,k3,indx,m); Ytemp[1] := Ytemp[1] + h; for i to m do Ytemp[i+1] := Ytemp[i+1] + R1 * k1[i] + R2 * k2[i] + R3 * k3[i] od; od; CC[laststored_cc] := CC[laststored_cc] + 1; for i from 0 to m do YY[CC[laststored_cc], i] := Ytemp[i + 1] od end: dointegration := proc(f, jac, YY, CC, Horg) local Echeck, X, DidX, DidY, dY, Eold, Converged, i; Echeck := 0; Eold := 0; Converged := false; DidX := array(1..round(CC[sMax_cc])); DidY := array(1..round(CC[m_cc]), 1..round(CC[sMax_cc])); dY := array(1..round(CC[m_cc])); while ((not Converged) and (CC[sUsed_cc] < CC[sMax_cc]) and (CC[Idid_cc] = 0) and (Echeck < CC[Emax_cc])) do if CC[sUsed_cc] > 0 then CC[h_cc] := Horg; CC[laststored_cc] := CC[laststored_cc] - 1; fi; CC[sUsed_cc] := CC[sUsed_cc] + 1; CC[h_cc] := Horg/`besirk/sSeq`[CC[sUsed_cc]]; if CC[method_cc] = 1 then `BESIRK/msirk3`(f, jac, YY, CC) elif CC[method_cc] = 2 then `BESIRK/rk4`(f, YY, CC) elif CC[method_cc] = 3 then `BESIRK/rk2`(f, YY, CC) fi; CC[ss_cc] := CC[ss_cc] + 1; X := CC[h_cc]; for i from 2 to CC[polord_cc] do X := X*CC[h_cc] od; `BESIRK/polpol`(CC, DidX, DidY, X, YY, dY); if CC[Idid_cc] = polerr then CC[Idid_cc] := 0; Echeck := CC[Emax_cc] fi; if ((CC[sUsed_cc] > 1) and (Echeck < CC[Emax_cc])) then Converged := StandT(CC, YY, dY); if ((not Converged) and (CC[sOld_cc] <> 0) and (CC[sUsed_cc] > CC[sOld_cc] + CC[sStep_cc])) then Echeck := CC[Emax_cc]; else if Eold = 0 then Eold := CC[Ecur_cc] else if (Eold > CC[Ecur_cc]) then Echeck := 0; Eold := CC[Ecur_cc]; else Echeck := Echeck + 1; fi; fi; fi; else Converged := false; fi; od; if Converged then YY[CC[laststored_cc], 0] := YY[CC[laststored_cc] - 1, 0] + Horg fi; if ((Echeck >= CC[Emax_cc]) or ((not Converged) and (CC[sUsed_cc] = CC[sMax_cc]))) then CC[Idid_cc] := fail; fi; if CC[Idid_cc] <> 0 then CC[laststored_cc] := CC[laststored_cc] - 1; CC[h_cc] := Horg; fi; end: findnexth := proc(CC, Horg) local sToUse, hnext, Shrink, Grow, MaxInc, i; Shrink := 0.9; Grow := 1.1; MaxInc := 10; if CC[Idid_cc] = 0 then sToUse := min(CC[sOld_cc] + CC[sStep_cc], CC[sUse_cc]); if (CC[sUsed_cc] = sToUse) then hnext := Shrink*Horg; elif (CC[sUsed_cc] = sToUse - 1) then hnext := Grow*Horg ; else if CC[sUsed_cc] = 0 then hnext := Horg*`besirk/sSeq`[sToUse - 1] else hnext := Horg*`besirk/sSeq`[sToUse - 1]/`besirk/sSeq`[CC[sUsed_cc]] fi; if hnext/Horg > MaxInc then hnext := Horg*MaxInc fi; fi; else if CC[Idid_cc] = fail then hnext := Horg/4.0; for i from 1 to round((CC[sMax_cc] - CC[sUse_cc])/2) do hnext := hnext/2.0 od; CC[Idid_cc] := 0 else hnext := 0; fi; fi; if hnext < Horg then CC[Reduct_cc] := 0 else if CC[Reduct_cc] = 0 then hnext := Horg fi; CC[Reduct_cc] := 1 fi; if CC[sOld_cc] = CC[sUsed_cc] then CC[sStep_cc] := 2 else CC[sStep_cc] := 1 fi; CC[sOld_cc] := CC[sUsed_cc]; if abs(hnext) < CC[hmin_cc] then CC[Idid_cc] := nexth fi; hnext; end: checkonh := proc(CC, t) local h; h := CC[h_cc]; if abs(h) > CC[hmax_cc] then h := CC[hmax_cc]*CC[Direc_cc] ; fi; if abs(h) > abs(CC[tout_cc] - t) then h := CC[Direc_cc]*abs(t - CC[tout_cc]); fi; if CC[hmin_cc] > 0 and abs(h) < CC[hmin_cc] then h := CC[hmin_cc]*CC[Direc_cc] ; fi; CC[h_cc] := h end: BESIRKshell := proc(f, jac, YY, CC) local hnext, iter, Horg; CC[Reduct_cc] := 1; CC[sStep_cc] := 1; CC[sOld_cc] := CC[sUse_cc]; iter := 0; CC[laststored_cc] := 0; while ((iter < CC[maxIt_cc]) and (CC[Idid_cc] = 0)) do iter := iter + 1; CC[sUsed_cc] := 0; CC[ss_cc] := 0; Horg := checkonh(CC, YY[CC[laststored_cc], 0]); dointegration(f, jac, YY, CC, Horg); hnext := findnexth(CC, Horg); if CC[res_cc] = 0 then if abs(YY[CC[laststored_cc], 0]) >= abs(CC[tout_cc]) then CC[Idid_cc] := ddone else if CC[Idid_cc] = 0 then CC[h_cc] := hnext fi fi; else if (abs(YY[CC[laststored_cc], 0]) >= abs(CC[tout_cc]) or YY[CC[laststored_cc],1] < CC[p_cc] or YY[CC[laststored_cc],1] > CC[q_cc] or YY[CC[laststored_cc],2] < CC[r_cc] or YY[CC[laststored_cc],2] > CC[s_cc]) then CC[Idid_cc] := ddone else if CC[Idid_cc] = 0 then CC[h_cc] := hnext; fi; fi; fi; od: end: besirkbandf := proc(f, jac, CCt, init) local YY, i, k, h, m, CC, Bfi, Bbi, precision, Bf, Bb, B, ti, tf; CC := array(1..48); CC := copy(CCt); ti := evalf(CC[ti_cc]); tf := evalf(CC[tf_cc]); m := round(CC[m_cc]); h := CC[h_cc]; YY := array(0..round(CC[maxIt_cc]), 0..m); precision := max(Digits, trunc(evalhf(Digits))); for i to m + 1 do YY[0, i - 1] := evalf(init[i], precision) od; Bfi := 0; Bbi := 0; if tf - init[1] > CC[hmin_cc] then CC[tout_cc] := tf; CC[Direc_cc] := 1; CC[h_cc] := h; if precision <= evalhf(Digits) then evalhf(BESIRKshell(f, jac, var(YY), var(CC))); else BESIRKshell(f, jac, YY, CC); fi; Bf := array(0..round(CC[laststored_cc])); Bfi := round(CC[laststored_cc]); for i from 0 to round(CC[laststored_cc]) do Bf[i] := [seq(YY[i, k], k = 0..m)] od; fi; CC := copy(CCt); if (evalf(init[1]) - evalf(ti)) > CC[hmin_cc] then CC[tout_cc] := ti; CC[Direc_cc] := -1; CC[h_cc] := -h; if precision <= evalhf(Digits) then evalhf(BESIRKshell(f, jac, var(YY), var(CC))); else BESIRKshell(f, jac, YY, CC); fi; Bb := array(0..round(CC[laststored_cc])); Bbi := round(CC[laststored_cc]); for i from 0 to round(CC[laststored_cc]) do Bb[i] := [seq(YY[i, k], k = 0..m)] od; fi; B := array(-Bbi .. Bfi); if Bbi > 0 then for i from -Bbi to 0 do B[i] := Bb[-i] od; fi; if Bfi > 0 then for i from 0 to Bfi do B[i] := Bf[i] od; fi; eval(B) end: besirk := proc(F) local B, CC, optargs, initvars, jac, m, f, init, func, vars; optargs := map( proc(xx) if type(xx,equation) then xx else 'dummy' = xx fi end, [args[2..nargs]]); optargs := ['intmethod' = 'besirk', 'numsteps' = 1, op(optargs)]; CC := array(1..48): initvars := `ODE/setup`(F, CC, optargs); vars := initvars[2]; init := initvars[1]; func := initvars[3]; f := `plots/BESIRK/makeproc`(func); if round(CC[method_cc]) = 1 then jac := `plots/ODE/makejac`(func) else jac := 1 fi; B := besirkbandf(f, jac, CC, init[1]); eval(B) end: orbitplot := proc(F) local B, CC, optargs, vv, ff, i, j, flows, backg, mystyle, jac, m, f, initlist, init, func, vars, initvars, ni, prange, p, q, r, s, u, v, pr, v1, v2, v3, nf, nb, viewo, prange2; optargs := map( proc(xx) if type(xx,equation) then xx else 'dummy' = xx fi end, [args[2..nargs]]); CC := array(1..48): initvars := `ODE/setup`(F, CC, optargs); vars := initvars[2]; initlist := initvars[1]; func := initvars[3]; vv := `ODE/graphics/setcc`(F, CC, optargs, vars); viewo := subs(optargs, 'view'); if viewo = 'view' then prange := [DEFAULT,DEFAULT] else if type(viewo[1],list) then viewo := viewo[1] fi; prange := [seq(subs(optargs,viewo[i]), i=1..nops(viewo))]; if nops(prange) = 2 then if not type(prange,[range, range]) then prange := [DEFAULT,DEFAULT] else p := op(1, op(1, prange)); q := op(2, op(1, prange)); r := op(1, op(2, prange)); s := op(2, op(2, prange)) fi elif nops(prange) = 3 then if not type(prange, [range, range, range]) then prange := [DEFAULT,DEFAULT, DEFAULT] else p := op(1, op(1, prange)); q := op(2, op(1, prange)); r := op(1, op(2, prange)); s := op(2, op(2, prange)); u := op(1, op(3, prange)); v := op(2, op(3, prange)) fi fi fi; prange2 := subs(optargs, 'plotrange'); if not prange2 = 'plotrange' then prange := prange2 fi; vv := vv[1]; f := `plots/BESIRK/makeproc`(func); if round(CC[method_cc]) = 1 then jac := `plots/ODE/makejac`(func) else jac := 1 fi; flows := NULL; i := 0; ni := nops(initlist); for init in initlist do i := i+1; if round(CC[intmethod_cc]) = 1 then B := besirkbandf(f, jac, CC, init); elif round(CC[intmethod_cc]) = 2 then B := `ODE/dsolve`(func, CC, init); else B := `ODE/rkhfbandf`(f, CC, init); fi; flows := flows, `ODE/makeflow`(B, CC, i, ni, vv, vars, optargs); od; backg := subs(select(xx -> (op(1,xx) = 'background'),optargs), background); if backg = 'background' then backg := NULL elif nops(prange) = 2 then backg := POLYGONS([[p,s],[p,r],[q,r],[q,s]], backg) else backg := NULL fi; pr := subs(select(xx -> (op(1,xx) = 'projections'),optargs), projections); if pr <> 'projections' then mystyle := `PATCHNOGRID`; pr := subs(pr, viewo); v1 := op(select(unapply(vars[xx] = viewo[1],xx) , [seq(j,j=1..nops(vars))])); v2 := op(select(unapply(vars[xx] = viewo[2],xx) , [seq(j,j=1..nops(vars))]));; v3 := op(select(unapply(vars[xx] = viewo[3],xx) , [seq(j,j=1..nops(vars))])); nf := op(2,op(2,eval(B))); nb := op(1,op(2,eval(B))); prange := [DEFAULT, DEFAULT, DEFAULT]; backg := POLYGONS([[pr[1],r,v],[pr[1],r,u],[pr[1],s,u],[pr[1],s,v]], COLOR(RGB,0.25,0.25,0.25)), POLYGONS([[p,pr[2],v],[p,pr[2],u],[q,pr[2],u],[q,pr[2],v]], COLOR(RGB,0.5,0.5,0.5)), POLYGONS([[q,r,pr[3]],[p,r,pr[3]],[p,s,pr[3]],[q,s,pr[3]]], COLOR(RGB,0.75,0.75,0.75)), CURVES([seq([pr[1],B[j][v2],B[j][v3]], j = nb .. nf)], black), CURVES([seq([B[j][v1],pr[2],B[j][v3]], j = nb .. nf)], black), CURVES([seq([B[j][v1],B[j][v2],pr[3]], j = nb .. nf)], black), ORIENTATION(-45,45) else mystyle := `HIDDEN` fi; if nops(prange) <> 3 then PLOT(VIEW(op(prange)), flows, backg, AXESSTYLE(BOX), STYLE(PATCHNOGRID)) else PLOT3D(VIEW(op(prange)), flows, backg, AXES(FRAME), STYLE(mystyle)) fi end: macro(_AA = _AA); macro(_body = _body); macro(_eq = _eq); macro(h_cc = h_cc); macro(m_cc = m_cc); macro(p_cc = p_cc); macro(q_cc = q_cc); macro(r_cc = r_cc); macro(s_cc = s_cc); macro(its_cc = its_cc); macro(steps_cc = steps_cc); macro(asp_cc = asp_cc); macro(fhead_cc = fhead_cc); macro(fishs_cc = fishs_cc); macro(fishl_cc = fishl_cc); macro(method_cc = method_cc); macro(sUse_cc = sUse_cc); macro(ss_cc = ss_cc); macro(hmin_cc = hmin_cc); macro(hmax_cc = hmax_cc); macro(Direc_cc = Direc_cc); macro(maxIt_cc = maxIt_cc); macro(Ecur_cc = Ecur_cc); macro(polord_cc = polord_cc); macro(laststored_cc = laststored_cc); macro(AbsTol_cc = AbsTol_cc); macro(RelTol_cc = RelTol_cc); macro(tout_cc = tout_cc); macro(sUsed_cc = sUsed_cc); macro(Emax_cc = Emax_cc); macro(sMax_cc = sMax_cc); macro(sOld_cc = sOld_cc); macro(sStep_cc = sStep_cc); macro(Reduct_cc = Reduct_cc); macro(a_cc = a_cc); macro(res_cc = res_cc); macro(numdiff_cc = numdiff_cc); macro(Idid_cc = Idid_cc); macro(hboxes_cc = hboxes_cc); macro(vboxes_cc = vboxes_cc); macro(intmethod_cc = intmethod_cc); macro(nclines_cc = nclines_cc); macro(flowf_cc = flowf_cc); macro(segments_cc = segments_cc); macro(orbthick_cc = orbthick_cc); macro(ti_cc = ti_cc); macro(tf_cc = tf_cc); macro(u_cc = u_cc); macro(v_cc = v_cc); macro(fthick_cc = fthick_cc); macro(field_cc = field_cc); macro(nreach = nreach); macro(fail = fail); macro(singul = singul); macro(userv = userv); macro(nexth = nexth); macro(ddone = ddone); macro(polerr = polerr); macro(node = node);