(*:Mathematica Version: 3.0 *) (*:Package Version: 1.0 *) (*:Name: DDSn` *) (*:Title: Exploring Discrete Dynamical Systems of Several Variable *) (*:Authors: Silvia Heubach and Jen Chen*) (*:Keywords: Iterative model equation, equilibrium, cobb-web graph, Eigenvalues and Eigenvectors *) (*:Summary: This package can be used in conjunction with the palette DDSnP. It works for discrete dynamical systems of the form x(n+1) = f(x(n)), where x(n) is the vector of the system variables. Functions compute system values from this equation, and the model values can be displayed as regular graph (with time on the horizontal axis) or as a graph of two variables. A animation of the meaning of Eigenvalues and Eigenvectors can be performed, and system dynamics can be visualized. The general solution for a linear DDS can be displayed for several cases *) BeginPackage["DDSn`","Graphics`Arrow`","Graphics`MultipleListPlot`"] IteratedValueSeqs::usage = "IteratedValueSeqs[fun,var,init,nmax,s] gives the values of the DDS defined by var(k+1)= fun(var(k)) for n = 0,...,nmax, in steps of size s, where var is the list of variables, fun is the list of functions defining the system, and init gives the value of the system at time 0. The step size s is optional; if no value is given, the step size will be 1. To get the value of the system only at time n = nmax, use step size s = nmax." SysEquil::usage = " SysEquil[{c1,c2,c3,c4}]computes the non-zero equilibrium point for the two-dimensional system x(n+1)= c1 x(n) + c2 x(n)y(n), y(n+1) = c3 y(n) + c4 x(n)y(n)" SysGraph::usage = "SysGraph[{c1,c2,c3,c4},init,nmax,joint,opts]displays the values of the two-dimensional system x(n+1)= c1 x(n) + c2 x(n)y(n), y(n+1) = c3 y(n) + c4 x(n)y(n) for n = 0,...,nmax, starting with the set of initial values given as pairs in the list init. If joint is True, then the values in each sequence are connected. The value(s) of opts indicate(s) ranges for the input and output values to be plotted given in the form xplot ->{xmin,xmax}, yplot -> {ymin, ymax}. These are optional, and can be used either separately or together." SysDyn::usage = "SysDyn[{c1,c2,c3,c4},range]displays the dynamics of the two-dimensional system x(n+1)= c1 x(n) + c2 x(n) y(n), y(n+1) = c3 y(n) + c4 x(n) y(n). If no range is given, only positive values of x and y are displayed. To see the dynamics for all possible values, use 'range -> All'." SeqsGraph::usage = "SeqsGraph[vals,vars,joint] displays the values of selected variables whose values were generated by the function IteratedValueSeqs and stored in the list vals. The list vars indicates which of the variables are to be displayed, where 0 = time, 1 = first variable, 2 = second variable,... If the optional value joint is set to True, the data points will be connected by line segments." LiveVecs::usage= "LiveVecs[A, num] creates a series of graphs of pairs of vectors {x, A.x}, and computes the Eigenvalues and Eigenvectors of the 2x2 matrix A. The value for num is optional and controls the number of different vectors x used; if no value is given, then num is set to 10." GenSol::usage= "GenSol[A] displays the general solution of the discrete linear system x(n+1) = A x(n), where x(n) is the vector of system variables" sol::usage= "sol[k] computes the general solution of the discrete linear system x(n+1) = A x(n) at time k, where x(k) is the vector of system variables." range::usage = "range determines the range of output values displayed in SysDyn" xplot::usage = "" Begin["`Private`"] Off[RuleDelayed::rhs]; Off[SetDelayed::write]; NewVals[fun_List,var_List,init_List, n_Integer] := NestList[Bigfunc[fun,var,#]&,init,n]; Bigfunc[fun_,var_,init_]:=N[fun /. Thread[var -> init]]; IteratedValueSeqs[fun_List,var_List,init_List,nmax_Integer,step_Integer:1]:= Module[{vals,seq,headline}, If[Length[fun]==Length[var]==Length[init], headline={"n"}; vals=NewVals[fun,var,init,nmax]; seq = Table[Flatten[{i-1,vals[[i]]}],{i,1,nmax+1,step}]; headline = Prepend[Table[DisplayForm [StringJoin[ToString[var[[i]]],"(n)"]], {i,1,Length[fun]}],"n"]; Prepend[seq,headline]//TableForm, Print["The number of functions must equal the number of variables and the number of initial values"]] ] (* End of Module *) ValSeqs[fun_List,var_List,init_List,nmax_Integer,step_Integer:1]:= Module[{vals,seq,headline}, headline={"n"}; vals=NewVals[fun,var,init,nmax]; seq = Table[Flatten[{i-1,vals[[i]]}],{i,1,nmax+1,step}]; headline = Prepend[Table[DisplayForm[StringJoin[ToString[var[[i]]],"(n)"]], {i,1,Length[fun]}],"n"]; Prepend[seq,headline]//TableForm ] (* End of Module *) SysEquil[{c1_,c2_,c3_,c4_}]:= Which[(c1 == 0 || c2 == 0 || c4 == 0 || c3 == 0), Print["All constants need to be non-zero"], c1 == 1, If[c3 == 1, Print["Both axes consist entirely of equilibrium points"],{"x",0}], c3 == 1,{0,"y"}, True, {(1-c3)/c4,(1-c1)/c2}] SysGraph[{c1_,c2_,c3_,c4_},init_,itsteps_, joint_, opts___?OptionQ]:= Module[{vals={},equilibrium, m = Length[init],vermin,vermax,hormin,hormax, xr,yr,start = init}, Clear[x,y]; h[{x_,y_}]={c1*x+c2*x*y,c3*y+c4*x*y}; If[(c1 == 0 || c2 == 0 || c4 == 0 || c3 == 0), Print["All constants need to be non-zero"], equilibrium={(1-c3)/c4,(1-c1)/c2}; Do[ AppendTo[vals,NestList[h,start[[i]],itsteps]],{i,m}]; (* Print[vals]; *) vermax = Max[Table[vals[[i,j,2]],{i,m},{j,itsteps+1}]]; vermin = Min[Min[Table[vals[[i,j,2]],{i,m},{j,itsteps+1}]],0]; hormax = Max[Table[vals[[i,j,1]],{i,m},{j,itsteps+1}]]; hormin = Min[Min[Table[vals[[i,j,1]],{i,m},{j,itsteps+1}]],0]; {xr,yr}= {xplot,yplot} /.{opts} /.{xplot -> {hormin, hormax},yplot -> {vermin,vermax}}; col ={RGBColor[1,0.7,0], RGBColor[0.3,0.8,0],RGBColor[1,0,1], RGBColor[0,0.9,1],GrayLevel[0.7]}; plotgraph = Table[ListPlot[vals[[i]], PlotJoined -> joint, PlotStyle -> {col[[Mod[i,4]+1]],PointSize[0.018]}, PlotRange -> {xr,yr},DisplayFunction->Identity],{i,1,m}]; graph5=Plot[equilibrium[[2]],{x,xr[[1]],xr[[2]]},PlotRange->yr,PlotStyle ->{RGBColor[1,0,0]}, DisplayFunction->Identity]; graph6=Show[Graphics[ { {RGBColor[1,0,0],Line[{{equilibrium[[1]],yr[[1]]},{equilibrium[[1]],yr[[2]]}}]}, {RGBColor[1,0,0],PointSize[0.02],Point[equilibrium]}, Table[ {col[[Mod[i,4]+1]],PointSize[0.02],Point[init[[i]]]},{i,1,m}], Table[Text[ToString[i],init[[i]]],{i,1,m}]}, DisplayFunction->Identity]]; Show[plotgraph,graph5,graph6, DisplayFunction->$DisplayFunction ]] ] (* End of Module *) Color={RGBColor[0,0,1] ,RGBColor[0.7,.3,0.5],RGBColor[1,0,0],RGBColor[0,1,0], RGBColor[1,0,1]}; SysDyn[{c1_,c2_,c3_,c4_},opts___?OptionQ]:= Module[ {xeq = (1-c3)/c4, yeq = (1-c1)/c2,minxeq,maxxeq,minyeq,maxyeq,xdist, ydist,plr}, minxeq = Min[0,xeq];maxxeq = Max[0,xeq]; minyeq = Min[0,yeq]; maxyeq = Max[0,yeq]; xdist = Abs[xeq];ydist = Abs[yeq]; xxarrow[{x_,y_}]:= If[c2>0, If[(x>0 && y >yeq) || (x < 0 && y < yeq),1,-1], If[(x>0 && y yeq),1,-1]]; yyarrow[{x_,y_}]:= If[c4>0, If[(y>0 && x >xeq) || (y< 0 && x< xeq),1,-1], If[(y>0 && x xeq),1,-1]]; arrowgraph[{x_,y_}]:= Graphics[{Arrow[{x,y},{x+xxarrow[{x,y}]*.25 xdist,y}, HeadLength->0.03], Arrow[{x,y},{x,y+yyarrow[{x,y}]*.25 ydist}, HeadLength->0.03]}]; If[(c1==1||c3==1), pts =Flatten[ Table[{i,j},{i,-1,1,2},{j,-1,1,2}],1], pts = Flatten[ Table[{i,j},{i,minxeq-.5 xdist,maxxeq+.5 xdist,xdist},{j, minyeq-.5ydist,maxyeq+.5ydist,ydist}],1]]; arrowlist = Map[arrowgraph,pts]; equilline = Graphics[{RGBColor[1,0,0],Line[{{minxeq-xdist,yeq},{maxxeq+xdist,yeq}}], Line[{{minxeq-xdist,0},{maxxeq+xdist,0}}], Line[{{xeq,minyeq-ydist},{xeq,maxyeq+ydist}}], Line[{{0,minyeq-ydist},{0,maxyeq+ydist}}],{PointSize[0.02], Point[{xeq,yeq}],Point[{0,0}]}}]; plr=range/.{opts}/.{range->{{0,maxxeq+xdist},{0,maxyeq+ydist}}}; Show[{arrowlist,equilline}, Axes ->True, PlotRange ->plr] ] (* End of Module *) SeqsGraph[name_,usernum_List,joint_:False]:=Module[ {ls, xlabel, ylabel, vertlabel, sel, gset, yorg, xorg, index}, sel = usernum; sortedsel = Sort[sel]; If[ListQ[name]==False,ls = name[[1]],ls = name]; If[(sortedsel[[1]]!=0 && Length[sel]>2)||(Max[sel] >= Length[Transpose[ls]]), Print["The list of variables to be graphed is incorrect"], (* else *) If[NumberQ[N[ls[[1,1]]]]==False, xlabel = ls[[1,sortedsel[[1]]+1]]; vertlabel = Table[ls[[1,sortedsel[[i]]+1]],{i,2,Length[sel]}]; ls = Transpose[Rest[ls]], (*else *) ls = Transpose[name]; If[sortedsel[[1]] == 0, xlabel = "time", xlabel = StringJoin["var",ToString[sortedsel[[1]]]]]; vertlabel = Table[StringJoin["var", ToString[sortedsel[[i]]]],{i,2,Length[sel]}]]; (*End If*) If[Length[sel]==2 && sortedsel[[1]]!=0, ylabel = vertlabel[[1]],ylabel = vertlabel]; If[sortedsel[[1]] == 0, index = 1, index = sortedsel[[1]]]; xorg = 1.1 Min[ls[[index]]] - .1 Max[ls[[index]]]; yorg = 1.1 Min[Map[Min, Rest[ls]]] - .1 Max[Map[Max,Rest[ls]]]; gset =Table[ ListPlot[Transpose[{ls[[sortedsel[[1]]+1]],ls[[sortedsel[[i]]+1]]}], PlotStyle-> {Color[[Mod[i,5]+1]], PointSize[0.02]},PlotJoined->joint, AxesLabel ->{xlabel,ylabel}, AxesOrigin -> {xorg,yorg}, DisplayFunction->Identity, PlotRange -> All], {i,2,Length[sel]}]; Show[{gset, Graphics[{Line[{{xorg,yorg},{ Min[ls[[index]]],yorg}}], Line[{{xorg,yorg},{xorg, Min[Map[Min, Rest[ls]]]}}]}]},DisplayFunction -> $DisplayFunction] ] (* End of If *) ] (* End of Module *) Norm[{x_,y_}]:=Sqrt[x^2+y^2] step[x_]:=1/2^(Floor[Log[2,10/x]]) LiveVecs[{{x1_,y1_},{x2_,y2_}},numframe_:15]:= Module[ {A={{x1,y1},{x2,y2}},lambda,eigenvecs,angles, diff = Floor[360/numframe],v1,v2,unitvecs,outvecs, maxlength,s,m,ringticks,rings}, (* I switch off the warning message when using the function ArcTan[] *) Off[ArcTan::indet]; Off[Part::partw]; lambda = N[Eigenvalues[A],3]; (* Error checking the Eigen values of the input matrix. If they are complex, then display an error message. *) If[Head[lambda[[1]]]===Complex || Head[lambda[[2]]]===Complex, Print[ "The Eigenvalues of your matrix are complex. Please choose \ another matrix which will produce real Eigenvalues."], (* else *) v1 =N[Eigenvectors[A][[1]]]; v2= N[Eigenvectors[A][[2]]]; angles = Table[i*diff,{i,0,numframe}]*Pi/180; If[v1[[2]]!=0, AppendTo[angles, Table[{ArcTan[v1[[1]],v1[[2]]], ArcTan[v1[[1]],v1[[2]]]+ Pi},{i,1,3}]], If[v1[[1]]!=0, AppendTo[angles, Table[{0,Pi},{i,3}]]]]; If[v2[[2]]!=0, AppendTo[angles, Table[{ArcTan[v2[[1]],v2[[2]]], ArcTan[v2[[1]],v2[[2]]]+Pi},{i,1,3}]], If[v2[[1]]!=0, AppendTo[angles, Table[{0,Pi},{i,3}]]]]; angles=Sort[N[Flatten[angles]]]; unitvecs = Chop[{Cos[angles],Sin[angles]}]; outvecs = A.unitvecs; maxlength = Max[Map[Norm,Transpose[outvecs]],1]; (* Now we get to the display *) s = step[maxlength]; m = Ceiling[maxlength]; ringticks = N[Complement[Range[-m,m,2*s],{0}],2]; rings = Graphics[{GrayLevel[.8],Table[Circle[{0,0},n],{n,s,m,s}]}]; InOutVecs[{{x_,y_},{u_,v_}}]:= Show[{rings, Graphics[{{RGBColor[1,0,0], Arrow[{0,0},{x,y},HeadLength -> 0.025]}, {RGBColor[0,0,1],Arrow[{0,0},{u,v},HeadLength -> 0.025]}}]}, AspectRatio -> Automatic, Axes ->True, Ticks ->{ringticks,ringticks}, PlotLabel -> TraditionalForm[StringForm["input = `1`", {x,y}]]] ; Print[DisplayForm[SubscriptBox["\[Lambda]","1"]],"= ",lambda[[1]]," ",DisplayForm[SubscriptBox["v","1"]]," = ",MatrixForm[v1/Norm[v1]]," " DisplayForm[SubscriptBox["\[Lambda]","2"]],"= ",lambda[[2]]," ",DisplayForm[SubscriptBox["v","2"]]," = ",MatrixForm[v2/Norm[v2]]]; Map[InOutVecs,Transpose[{Transpose[unitvecs],Transpose[outvecs]}]] ] (* End If *) ] (* End Module *) DoubleRootDep[mat_?MatrixQ]:= Module[{values,sys,lambda,vec1,col={x,y},num1,num2,secondvec}, sys=Eigensystem[mat]; lambda=sys[[1,1]]; If[sys[[2,1]]=={0,0},vec1=sys[[2,2]],vec1=sys[[2,1]]]; values=Solve[(mat-lambda*IdentityMatrix[2]).col==lambda*vec1]; If[Length[values]==2,num1=values[[1,1,2]]; num2=values[[1,2,2]], If[values[[1,1,1]]===x, num1=values[[1,1,2]];num2=1, num1=1; num2=values[[1,1,2]]] (* End of IF *) ]; (* End of IF *) secondvec={num1,num2}; Print["x(k) = ", DisplayForm[SubscriptBox["c","1"]], "*",DisplayForm[SuperscriptBox[lambda,"k"]], "*",MatrixForm[vec1]," + ", DisplayForm[SubscriptBox["c","2"]], "*{k*", DisplayForm[SuperscriptBox[lambda,"k"]], "*",MatrixForm[vec1]," + ", DisplayForm[SuperscriptBox[lambda,"k"]],"*", MatrixForm[secondvec],"}" ]; (* End of Print *) sol[k_]:=Subscript["c",1]*lambda^k*vec1 + Subscript["c",2]*(k*lambda^k*vec1 + lambda^k*secondvec) ] (* End of Module *) GenSol[mat_?MatrixQ]:= Module[{Eval,Evec,num,mat1=mat,result,realpos,impos,imvals,imvecs,realvals, realvecs,str,solterms,newstr,foo,angle,lamb1,lamb2,vec1,vec2,rad}, Clear[sol]; num=Length[mat]; If[num<2, Print["You do not have a problem with 2 or more variables"], (* else *) Eval=Eigenvalues[mat]; Evec= Eigenvectors[mat]; (* Checking for dependence of EigenVectors *) If[Length[NullSpace[Transpose[Evec]]] > 0, If[num > 2,Print["Can not compute general solution"], (* Double root with only one eigenvector *) DoubleRootDep[mat1]], (* Vectors are all independent *) realpos = Position[Im[Eval],0]; impos=Partition[ Complement[Range[num],Flatten[realpos]],1]; (* Print["realpos = ",realpos," impos = ",impos]; *) If[impos=={}, (* no complex value *) realvals=Eval; realvecs=Evec, (* else *) realvals=Delete[Eval,impos]; realvecs=Delete[Evec,impos]; imvals=Delete[Eval,realpos]; imvecs=Delete[Evec,realpos]; (*Print["realvals = ",realvals, " imvals = ",imvals,"\n realvecs = ", realvecs, " imvecs = ",imvecs];*) (* Solving Complex cases *) For[l=1, l< Length[imvals], l=l+2, If[Re[imvals[[l]]] == Re[imvals[[l+1]]] && Im[imvals[[l]]] == -Im[imvals[[l+1]]], lamb1 = Re[imvals[[l]]]; lamb2 = Abs[Im[imvals[[l]]]]; vec1 = Simplify[(imvecs[[l]]+imvecs[[l+1]])/2]; vec2 = Simplify[(imvecs[[l]]-imvecs[[l+1]])/(2I)]; angle = N[ArcTan[lamb1,lamb2]]; rad = Sqrt[lamb1^2+lamb2^2]; AppendTo[realvals,{rad,rad}]; AppendTo[realvecs, (vec1*Cos[k*angle]+vec2*Sin[k*angle])]; AppendTo[realvecs, (vec2*Cos[k*angle]- vec1*Sin[k*angle])], Print["Something is wrong with Mathematica's ordering"] ]; (* End of IF *) ] (* End of FOR *) ]; (* End of IF *) realvecs = realvecs /. {TestPac`Private`k -> Global`k}; realvals=Flatten[realvals]; (*Print["realvals = ",realvals," realvecs = ",realvecs];*) Clear[str]; str={"x(k) = "}; solterms=Table[{SubscriptBox["c",ToString[i]],"*", SuperscriptBox[ToBoxes[realvals[[i]]],"k"], "*",ToBoxes[Partition[realvecs[[i]],1], TraditionalForm],"+"},{i,1,num}]; newstr=Drop[Flatten[AppendTo[str,solterms]],-1]; foo=RowBox[newstr] // DisplayForm; Print[foo]; cons=Table[Subscript["c",i],{i,1,num}]; sol[x_]:=(Apply[Plus,cons*realvals^x*realvecs]/. {Global`k -> x}) ] ] ] (* End of Module *) End[] Protect[LiveVecs,ItVals,ValSeqs,GenSol] EndPackage[]