{VERSION 5 0 "IBM INTEL NT" "5.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Headi ng 1" -1 3 1 {CSTYLE "" -1 -1 "Times" 1 18 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 4 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 2" -1 4 1 {CSTYLE "" -1 -1 "Times" 1 14 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 2 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 3" -1 5 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 1 1 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Author" 0 19 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 8 8 0 0 0 0 0 0 -1 0 }{PSTYLE "Title" -1 256 1 {CSTYLE "" -1 -1 "Times " 1 24 0 0 0 1 2 1 1 2 2 2 1 1 1 1 }3 1 0 0 12 12 1 0 1 0 2 2 19 1 }} {SECT 0 {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 256 "" 0 "" {TEXT -1 39 "Integrating the Frenet-Serret equations " }}{PARA 19 "" 0 "" {TEXT -1 168 "Modified from Matthias Kawski \+ http://math.asu.edu/~kawski\nDept. of Mathematics and Statis tics kawski@asu.edu\nArizona State University " }}} {SECT 1 {PARA 3 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 15 "Author and Date" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 "\n " }}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 24 "Content, Purpose and Use" } }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 0 {PARA 4 "" 0 " " {TEXT -1 32 "Updates and log of modifications" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 205 "Completely rewritten on Feb 3, 2000\n\nJanuary 2003 : Updated to MAPLE V release 8 -- only checked that it is running \n \+ as expected, but no effort to update e.g. linalg \+ to LinearAlgebra." }}}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 48 "Constant \+ curvature and torsion -- using linalg()" }}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 49 "Create the matrix A and calculate its exponential" }} {SECT 0 {PARA 5 "" 0 "" {TEXT -1 30 "Create the matrix A(kappa,tau)" } }{EXCHG {PARA 0 "" 0 "" {TEXT -1 140 "In the case of a constant coeffi cients there is no need to resort to using\ngeneral differential equat ion solvers -- linear algebra suffices." }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "restart;\nwith(linalg):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "assume(kappa,real);\nassume(tau,rea l);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 79 "A:=matrix(3,3,0); \nA[1,2]:=-kappa;\nA[2,1]:= kappa;\nA[2,3]:=-tau;\nA[3,2]:= tau;\n\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "print(A);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 0 {PARA 5 "" 0 "" {TEXT -1 39 "Optional detour 1: Via eigenvectors...." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "ev:=eigenvects(A);\n" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 35 "DD:=diag(seq(op(1,ev[k]),k=1..3));\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "Uinv:=blockmatrix(1,3,[seq(op(op(3, ev[k])),k=1..3)]);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "sim plify(multiply(htranspose(Uinv),Uinv));\n" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 53 "U:=inverse(Uinv):\nmap(simplify,multiply(Uinv,DD,U) );\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "expDD:=map(q->exp(s *q),DD);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "expsA:=map(si mplify,multiply(Uinv,expDD,U));\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 0 {PARA 5 "" 0 "" {TEXT -1 49 "Optional de tour 2: Using the Jordan decomposition" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 243 "The command jordan() basically calculates the same as ei genvects, but returns\nthe results already assembled in matrix format. Subtile differences may involve\nthe handling of repeated (and/or com plex) eigenvalues and rank deficient eigenspaces." }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "DD:=jordan(A,'Uinv');" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "print(Uinv);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "U:=inverse(Uinv):\nmap(simplify,mul tiply(Uinv,DD,U));\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "\ne xpsA:=map(simplify,multiply(Uinv,map(q->exp(s*q),DD),U));\n" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 0 {PARA 5 "" 0 " " {TEXT -1 40 "The standard route: Use exponential(...)" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 235 "In general it is safer -- fewer errors e tc. -- to simply let MAPLE do all the\nsmall steps presented in the tw o prior subsections. The command eponential(...)\ncalculates the matr ix exponential (fundamental matrix solution) in one step!" }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "'A'=evalm(A);\n " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "expsA:=exponential(A,s) :\nexp(s*'A')=evalm(expsA);\n\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 60 "Integrate the compl ete initial value problem -- general form" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "R0:=matrix(3,3,(i,j)->R||i||j||0);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "gamma0:=[x0,y0,z0];\n" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 78 "then the (T,N,B) frame at time s equals - - we only display one matrix entry --" }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "R:=multiply(R0,expsA):\nR[1,1];" }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 252 "and the curve requires one addit ional integration -- note, s is used as dummy integration\nvariable, a nd hence gamma is written as an expression depending on t rather than \+ on the \nmore natural parameter s.... -- again we only dispaly the fir st component" }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "gammas:=[seq(gamma0[k]+int(R[k,1],s=0..t),k=1..3)]:\ngammas[1] ;" }}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 18 "A specific example" }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "helixgeom := kappa = 2/5, ta u = 1/5;\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "expsA1:=map(q ->convert(subs(helixgeom,q),trig),expsA);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "helixIV:=x0=2,y0=0,z0=0;\n" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 83 "R01:=transpose(matrix(3,3,[0,2/sqrt(5),1/sqrt( 5),-1,0,0,0,-1/sqrt(5),2/sqrt(5)]));\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "R1:=multiply(R01,expsA1);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "gamma1:=[seq(subs(helixIV,gamma0[k])+int(R1[k,1] ,s=0..t),k=1..3)];\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 868 "wi th(plots):\nTfin:=4*Pi*sqrt(5):\ncurvepix:=spacecurve(gamma1,t=0..Tfin ,\n axes=normal,labels=[`x`,`y`,`z`],\n thic kness=3,shading=ZHUE,\n tickmarks=[2,2,2],orientation=[40 ,60]):\nNframes:=24:\ndisplay([seq(display([curvepix,\n \+ spacecurve([gamma1,gamma1+3*subs(s=t,[seq(R1[j,1],j=1..3)])],\n \+ color=blue,thickness=3),\n spacec urve([gamma1,gamma1+3*subs(s=t,[seq(R1[j,2],j=1..3)])],\n \+ color=green,thickness=3),\n spacecurve([ gamma1,gamma1+3*subs(s=t,[seq(R1[j,3],j=1..3)])],\n \+ color=red,thickness=3) \n ]),\n \+ t=[seq(k*Tfin/Nframes,k=0..Nframes)])],\n insequence=true, \n title=`Observe that N is always horizontal`,\n scalin g=constrained,orientation=[45,75]);\n\nt:='t':\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 48 "Variable curvatur e and torsion -- using dsolve()" }}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 55 "Create the general data structures to be used by dsolve" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 76 "Create the set of differential equations \+ -- here without using any matrices." }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "for i to 3 do\n for j to 3 do\n \+ a||i||j(s):=0;\n od;\nod;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 79 "a12(s) := -kappa(s);\na21(s) := kappa(s);\na23(s) := -tau(s);\na32(s) := tau(s);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 141 "DEs:=\{\n seq(seq(\n diff(R||i||j(s),s)=\n \+ convert([seq(R||i||k(s)*a||k||j(s),k=1..3)],`+`),\n j=1..3),i= 1..3)\n \};" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 127 "delta:= (i,j)->if (i=j) then 1 else 0 fi:\nICs:=\{\n seq(seq(\n R| |i||j(0)=delta(i,j),\n j=1..3),i=1..3)\n \};\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 76 "vars:=\{\n seq(seq(\n R ||i||j(s),\n j=1..3),i=1..3)\n \};\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 54 "Tes t: Again integrate the Frenet system for the helix " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "curve1:=kappa(s)=2/5,tau(s)=1/5;\n" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "DEs1:=subs(curve1,DEs);\n" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "sol:=dsolve(DEs1 union ICs , vars);\n" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 227 "This appears to be OK -- as expected, dsolve can handle linear constant coefficient\nsys tems -- but it took a littlelonger than just using linear algebra \"by hand\".\nFor convenience we may write the solution again in matrix fo rm:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 75 "with(linalg):\nR1:=m atrix(subs(sol,[seq([seq(R||i||j(s),j=1..3)],i=1..3)]));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 84 "TNB0:=transpose(matrix(3,3,[0,2/sqr t(5),1/sqrt(5),-1,0,0,0,-1/sqrt(5),2/sqrt(5)]));\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "TNB:=multiply(TNB0,R1);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "gamma0:=[2,0,0]:\ngamma1:=[seq(gamma0[k]+ int(TNB[k,1],s=0..t),k=1..3)];\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 818 "with(plots):\nTfin:=4*Pi*sqrt(5):\ncurvepix:=spacecu rve(gamma1,t=0..Tfin,\n axes=normal,labels=[`x`,`y`,`z`], \n thickness=3,shading=ZHUE,\n tickmarks=[2, 2,2],orientation=[40,60]):\nNframes:=24:\ndisplay([seq(display([curvep ix,\n spacecurve([gamma1,gamma1+3*subs(s=t,[seq(TNB[j,1 ],j=1..3)])],\n color=blue,thickness=3),\n \+ spacecurve([gamma1,gamma1+3*subs(s=t,[seq(TNB[j,2],j=1. .3)])],\n color=green,thickness=3),\n \+ spacecurve([gamma1,gamma1+3*subs(s=t,[seq(TNB[j,3],j=1..3)]) ],\n color=red,thickness=3) \n \+ ]),\n t=[seq(k*Tfin/Nframes,k=0..Nframes)])],\n \+ insequence=true,\n scaling=constrained,orientation=[45,75 ]);\n\nt:='t':\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}} {SECT 0 {PARA 4 "" 0 "" {TEXT -1 60 "A first variable curvature exampl e -- and a numeric solution" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 118 "As expected, dsolve cannot find a general solution in nice symbolic form . Thus\nrevert to ask for a numerical solution." }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 47 "sol:=dsolve(DEs1 union ICs, vars,type=numeric) ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "sol(0);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 95 "with(linalg):\nMsol:=t->matrix(subs (op(sol(t)[2..10]),[seq([seq(R||i||j(s),j=1..3)],i=1..3)])):\n" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 93 "Msol(2);\nTNB0:=transpose(ma trix(3,3,[0,2/sqrt(5),1/sqrt(5),-1,0,0,0,-1/sqrt(5),2/sqrt(5)]));\n" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "TNB:=t->map(evalf,multiply (TNB0,Msol(t))):\nTNB(2);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "ortho:=t->multiply(transpose(TNB(t)),TNB(t)):\northo(10);\n" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 288 "distort:=t->map(op,convert( matadd(ortho(t),diag(-1,-1,-1)),listlist)):\ndistort(10);\ndistort(100 );\ndistorts:=[seq(distort(2*j),j=0..20)]:\nwith(plots):\ndisplay([seq (plot(\n [seq([2*j-2,distorts[j,k]],j=1..21)],\n thick ness=2),\n k=1..9)],title=`Increasing errors`);\n" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 174 "Rather than again using dsolve wi th the above numerical solution as input to find \nthe associated curv e starts from the beginning, and integrates all 12 Des at the\nsame ti me." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 46 "The enlarged set of DEs -- all twelve together" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 76 "C reate the set of differential equations -- here without using any matr ices." }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 183 "for i to 3 do\n for j to 3 do\n a||i||j(t):=0;\n od; \nod:\na12(t) := -speed(t)*kappa(t);\na21(t) := speed(t)*kappa(t);\na 23(t) := -speed(t)*tau(t);\na32(t) := speed(t)*tau(t);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 199 "DEs:=\{\n seq(seq(\n d iff(R||i||j(t),t)=\n convert([seq(R||i||k(t)*a||k||j(t),k=1..3 )],`+`),\n j=1..3),i=1..3),\n seq(diff(gamma||k(t),t)=spee d(t)*R||k||1(t),k=1..3)\n \};" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 159 "delta:=(i,j)->if (i=j) then 1 else 0 fi:\nICs:=\{\n \+ seq(seq(\n R||i||j(0)=delta(i,j),\n j=1..3),i=1..3 ),\n seq(gamma||k(0)=0,k=1..3)\n \};\n" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 106 "vars:=\{\n seq(seq(\n R||i||j(t), \n j=1..3),i=1..3),\n seq(gamma||k(t),k=1..3)\n \};\n " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 0 {PARA 4 " " 0 "" {TEXT -1 29 "Creating nontrivial test data" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 296 "viviani:=[(1+cos(2*t)),sin(2*t),2*sin(t)]:\n with(plots):\ndisplay([sphereplot(1.95,0..2*Pi,0..Pi,style=wireframe,s hading=ZHUE),\n spacecurve(viviani,t=0..2*Pi,thickness=4,color =black)],\n scaling=constrained,axes=normal,tickmarks=[0,0,0],\n \+ orientation=[55,77],title=`Viviani's curve`);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 112 "v:=diff(viviani,t);\n#v:=diff([2*cos(t/sqrt( 5)),2*sin(t/sqrt(5)),t/sqrt(5)],t);\na:=diff(v,t);\naprime:=diff(a,t); \n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 210 "with(linalg):\nsp:=s implify(sqrt(dotprod(v,v,orthogonal)));\nBB:=simplify(crossprod(v,a)); \nbb:=simplify(sqrt(dotprod(BB,BB,orthogonal)));\ncur:=simplify(bb/sp^ 3);\ntor:=simplify(dotprod(BB,aprime,orthogonal)/bb^2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 163 "start:=eval(viviani,t=0);\nTini:=s implify(eval(map(q->q/sp,v),t=0));\nBini:=simplify(eval(map(q->q/bb,BB ),t=0));\nNini:=simplify(convert(crossprod(Bini,Tini),list));\n" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 0 {PARA 4 "" 0 " " {TEXT -1 45 "Numerical solution of all twelve DEs, example" }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 143 "curve2:=speed(t)=2*sqrt(1+(cos(t))^2),\n ka ppa(t)=sqrt(5+3*(cos(t))^2)/2/(1+cos((t))^2)^(3/2),\n tau(t)=3* cos(t)/(5+3*(cos(t))^2);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 97 "plot(subs(curve2,[speed(t),kappa(t),tau(t)]),t=0..2*Pi,\n colo r=[black,red,blue],thickness=3);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "DEs2:=subs(curve2,DEs);\n" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 219 "ICs2:=\{R11(0)=0, R12(0)=-1,R13(0)=0,\n \+ R21(0)=sqrt(2)/2, R22(0)= 0,R23(0)=-sqrt(2)/2,\n R31(0)=sqrt (2)/2, R32(0)= 0,R33(0)= sqrt(2)/2,\n gamma||1(0)=2,\n gam ma||2(0)=0,\n gamma||3(0)=0\};" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "sol:=dsolve(DEs2 union ICs2, vars,type=numeric);\n" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "sol(2);\n" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 136 "Outputs from dsolve(...,numeric) are usually p lotted with odeplot -- but here we do this\nonly for the curve, not fo r the Frenet frame..." }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 397 "display([sphereplot(1.95,0..2*Pi,0..Pi,style=wirefra me,\n grid=[24,18],shading=ZHUE),\n spacecu rve(viviani,t=0..2*Pi,thickness=2,color=blue),\n odeplot(sol,[ seq(gamma||k(t),k=1..3)],0..2*Pi,\n thickness=3,color=b lack,numpoints=25)],\n scaling=constrained,axes=normal,tickmarks=[0 ,0,0],\n orientation=[55,77],title=`Recreation of Viviani's curve`) ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 268 "display([plot([seq([t ,viviani[k],t=0..2*Pi],k=1..3)],\n color=[black,blue,mage nta],thickness=1),\n odeplot(sol,[seq([t,gamma||k(t)],k=1..3)] ,0..2*Pi,\n thickness=3,numpoints=25)\n ],title =`Comparison of curve w/ original data`);" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 55 "[seq(subs(sol(tk),[seq(R||k||3(t),k=1..3)]),tk=0..1 )];\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1171 "sph:=sphereplot( 1.95,0..2*Pi,0..Pi,style=wireframe,\n grid=[24,18], color=cyan):\ncurvepix:=odeplot(sol,[seq(gamma||k(t),k=1..3)],0..2*Pi, \n thickness=3,color=black,numpoints=100):\nNframes:=48 :\ndisplay([seq(display([curvepix,sph,\n spacecurve([[s eq(subs(sol(tk),gamma||k(t) ),k=1..3)],\n \+ [seq(subs(sol(tk),gamma||k(t)+R||k||1(t)),k=1..3)]],\n \+ color=blue,thickness=3),\n spacecurv e([[seq(subs(sol(tk),gamma||k(t) ),k=1..3)],\n \+ [seq(subs(sol(tk),gamma||k(t)+R||k||2(t)),k=1..3)]],\n \+ color=red,thickness=3),\n space curve([[seq(subs(sol(tk),gamma||k(t) ),k=1..3)],\n \+ [seq(subs(sol(tk),gamma||k(t)+R||k||3(t)),k=1..3)]],\n color=green,thickness=3)\n \+ ]),\n tk=[seq(k*2*Pi/Nframes,k=0..Nframes)])],\n \+ insequence=true,\n axes=normal,labels=[`x`,`y`,`z`],tic kmarks=[2,2,2],\n scaling=constrained,orientation=[55,75],\n \+ title=`Integrating speed, curv, and torsion`);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 25 "One proced ure does it all" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 43 "Define the procedure (execute \+ this section)" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2693 "intcurv:= proc()\n\nlocal speed, kappa, tau, trange, DEs, ICs, RR, xyz,\n i , j, k, curvepix, vars, sol, Nframes, mytitle;\nwith(plots):\nif nargs > 0 \n then kappa:=args[1]\n else \n kappa:=sqrt(5+3*(cos(t))^2 )/2/(1+cos((t))^2)^(3/2);\n fi;\nif nargs > 1 \n then tau:=args[2] \n else \n tau:=3*cos(t)/(5+3*(cos(t))^2);\n fi;\nif nargs > 2 \+ \n then speed:=args[3]\n else \n speed:=2*sqrt(1+(cos(t))^2);\n \+ fi;\nif nargs > 3 \n then trange:=args[4]\n else \n trange:=0. .2*Pi;\n fi;\nif nargs > 4 \n then \n RR:=args[5]\n else RR:=[ [0,sqrt(2)/2,sqrt(2)/2],\n [-1,0,0],\n [0,-sqr t(2)/2,sqrt(2)/2]]\n fi;\nif nargs > 5 \n then \n xyz:=args[6]\n else \n xyz:=[2,0,0] \n fi;\nif nargs > 6 \n then \n Nframe s:=args[7]\n else Nframes:=48;\n fi;\nif nargs > 7 \n then \n \+ mytitle:=args[8]\n else mytitle:=` `;\n fi;\nICs:=\{R11(0)=RR[1,1] , R12(0)= RR[2,1],R13(0)= RR[3,1],\n R21(0)=RR[1,2], R22(0)= RR[2 ,2],R23(0)= RR[3,2],\n R31(0)=RR[1,3], R32(0)= RR[2,3],R33(0)= RR [3,3],\n gamma||1(0)=xyz[1],\n gamma||2(0)=xyz[2],\n \+ gamma||3(0)=xyz[3]\};\n\nfor i to 3 do\n for j to 3 do\n a| |i||j(t):=0;\n od;\nod:\na12(t) := -speed*kappa;\na21(t) := speed* kappa;\na23(t) := -speed*tau;\na32(t) := speed*tau;\nDEs:=\{\n se q(seq(\n diff(R||i||j(t),t)=\n convert([seq(R||i||k(t) *a||k||j(t),k=1..3)],`+`),\n j=1..3),i=1..3),\n seq(diff(g amma||k(t),t)=speed*R||k||1(t),k=1..3)\n \};\nvars:=\{\n seq(s eq(\n R||i||j(t),\n j=1..3),i=1..3),\n seq(gamma|| k(t),k=1..3)\n \};\nsol:=dsolve(DEs union ICs, vars,type=numeric); \ncurvepix:=odeplot(sol,[seq(gamma||k(t),k=1..3)],trange,\n \+ thickness=3,color=black,numpoints=100):\ndisplay([seq(display([cu rvepix,\n spacecurve([[seq(subs(sol(tk),gamma||k(t) \+ ),k=1..3)],\n [seq(subs(sol(tk),gamma| |k(t)+R||k||1(t)),k=1..3)]],\n color=blue,t hickness=3),\n spacecurve([[seq(subs(sol(tk),gamma||k(t ) ),k=1..3)],\n [seq(subs(sol(tk),g amma||k(t)+R||k||2(t)),k=1..3)]],\n color=r ed,thickness=3),\n spacecurve([[seq(subs(sol(tk),gamma| |k(t) ),k=1..3)],\n [seq(subs(sol(t k),gamma||k(t)+R||k||3(t)),k=1..3)]],\n col or=green,thickness=3)\n ]),\n tk=[seq (op(1,trange)+k*\n (op(2,trange)-op(1,trange))/Nfra mes,k=0..Nframes)])],\n insequence=true,\n axes=norm al,labels=[`x`,`y`,`z`],tickmarks=[2,2,2],\n scaling=constrai ned,orientation=[55,75],\n title=mytitle);\nend:\n" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 0 {PARA 4 "" 0 " " {TEXT -1 8 "Examples" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "in tcurv();\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "intcurv(0,0,1 ,0..3);\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "intcurv(2/3,0) ;\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 154 "intcurv(2/5,1/5,1,0 ..8*Pi,\n [[0 ,2/sqrt(5), 1/sqrt(5)],\n [-1,0, \+ 0 ],\n [0 ,-1/sqrt(5),2/sqrt(5)]],\n [2,0,0]);\n " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 128 "intcurv(t*sin(t)/3,0,1 ,-3*Pi..3*Pi,\n [[1,0,0],[0,1,0],[0,0,1]],\n [0,0,0],8 0,\n `View me from the top!`);\n" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 12 "Fun in space" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 88 "intcurv(t*sin(t),1,1,-3*Pi..3*Pi,\n [[1,0,0],[0,1,0],[0,0, 1]],\n [0,0,0]);" }}}}}}{MARK "1 1 0" 168 }{VIEWOPTS 1 1 0 3 2 1804 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }