## dharr

Dr. David Harrington

## 6337 Reputation

19 years, 352 days
University of Victoria
Professor or university staff

## Social Networks and Content at Maplesoft.com

I am a professor of chemistry at the University of Victoria, BC, Canada, where my research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

## from eqns to Matrix...

I made the opposite assumption to @acer (that you started with the equations), and kept implicitplot butdid them all in one, which is simple here despite the inefficiency.

 >

 >

Define the equations to solve and plot

 >

Convert to the corresponding Matrix of coefficients and Vector (because we shouldn't have to enter equivalent information twice, and less likely to make a mistake)

 >

 >

We can do all equations in one implicitplot. (I like the brighter colors, so I left the quotes off.)

 >

 >

## Fourier-Bessel...

Depends a bit on how you want to enter it.

 >
 > A := m -> int(xi*(1 - xi^2)*BesselJ(0, BesselJZeros(0, m)*xi), xi = 0 .. 1, numeric)/           int(xi*BesselJ(0, BesselJZeros(0, m)*xi)^2, xi = 0 .. 1, numeric);

 >

 >

 >

Save an integral - see DLMF

 > A := m -> 2*int(xi*(1 - xi^2)*BesselJ(0, BesselJZeros(0, m)*xi), xi = 0 .. 1, numeric)/           evalf(BesselJ(1, BesselJZeros(0, m))^2);

 >

 >

## BesselJZeros...

`evalf(BesselJZeros(0, 1..5));`

gives: 2.404825558, 5.520078110, 8.653727913, 11.79153444, 14.93091771

## Minpoly...

Here's how I would do it. Not sure why I get the two different answers.

 >
 >

 >

 > use alpha=1-(1/2)/(1-(RootOf(16*_Z*(_Z*(2*_Z*(_Z*(8*_Z*(_Z*(_Z*(_Z*(32*_Z*(8*_Z-33)+1513)-812)-13)+267)-1469)-330)+811)+279)+345,index=2)-1/2)**2) in         expr:=(1+alpha)*sqrt(1-alpha**2)+(3+4*alpha)/12*sqrt(3-4*alpha**2)+2*(1+alpha)/3*sqrt(2*(1+alpha)*(1-2*alpha))+(1+2*alpha)/6*sqrt(2*((1-alpha)**2-3*alpha**2)) end:
 >

Looks correct

 >

 >

 >

Square roots problematic, so convert to RootOf

 >

We want the minimum polynomial. Over the rationals - not sure why this is not the same as above

 >

With as a field extension it is a quadratic

 >

 >

## Matrix assignments...

You can assign part of a Matrix to part of another one using the notation for row and column selection, for example

`A[6..7, 8..11] := M[3..4, 6..9]`

copies rows 3..4 and colums 6..9 of M into rows 5..7, cols 8..11 of A. More complicared variations are possible, in which the rows and columns do not have to be adjacent - see the MVextract help page

Edit: I changed the above so the blocks are the same size. Orignally I had A[5..7, 8..11] := M[3..4, 6..9] - see Joe Riel for what happens in this case.

## unfolding...

Here's a variation of @mmcdara's idea for a tetrahedron, since he challenged me :). Would require significantly more work for shapes where the folds aren't all on a flat base.

 > restart;
 > with(geom3d): with(GraphTheory):

Example - tetrahedron

 > solid := tetrahedron(t, point(o, [0,0,0]), 3); f:=faces(solid); pts:=vertices(solid); seq(point(cat(p,i),pts[i]),i=1..nops(pts));

Code the faces with the indices of the points.

 > indextable:=table(pts=~[\$1..nops(pts)]): ff:=subsindets(f,[algebraic\$3],x->indextable[x]);

Make corresponding graph

 > triangles:=map(CycleGraph,ff): tetgraph:=Graph([\$1..nops(pts)],`union`(map(Edges,triangles)[])): SetVertexPositions(tetgraph,evalf(pts)):

Try some cuts (red edges on drawing) - spanning tree works for the convex platonic solids

 > tree:=SpanningTree(tetgraph): HighlightSubgraph(tetgraph,tree);
 > DrawGraph(tetgraph, dimension = 3, size = [300, 300],orientation=[50,70]); #red edges are cuts

Uncut edges are the folds, which in this case are all in the same plane.

 > folds:=Edges(tetgraph) minus Edges(tree); basevertices:=convert(`union`(%[]),list); plane(base,map2(cat,p,basevertices)): triangle(basetriangle,map2(cat,p,basevertices)):

Remove base from the other faces

 > tofold:=remove(x->{x[]}={basevertices[]},ff);

Unfold

 > for m to nops(folds) do  fold:=folds[m];  foldface:=select(x-> fold subset {x[]}, tofold)[]; # find face to fold  facepts:=map2(cat,p,foldface); # points in face  angle:=evalf(FindAngle(base,plane(pl,facepts))); # angle between base and face  ii := select(has,foldface,fold); # fold sorted as in tofold  line(cat(raxis,m),map2(cat,p,ii)); # fold is rotation axis  rotation(cat(unfolded,m),triangle(cat(tr,m),facepts),Pi-angle,cat(raxis,m)); # and rotate end do:

Unfolded

 > draw([base,unfolded1,unfolded2,unfolded3]);

Angles all the same so generate animation

 > n:=20: angles:=[seq(0..Pi-angle,numelems=n)]:
 > for i to nops(angles) do   ang:=angles[i];   rotation(rot1,tr1,ang,raxis1);   rotation(rot2,tr2,ang,raxis2);   rotation(rot3,tr3,ang,raxis3);   dr[i]:=draw([solid,rot1,rot2,rot3]); end do:
 > plots:-display([seq(dr[i],i=1..n)],insequence=true)

 >

## matrix way...

Here's a more matrixy way.

```MM:=Matrix(2,[0,1,1,0]);
eq14 := lhs(eq13) = MM% . simplify(eval(MM^(-1) . rhs(eq13), `%.` = `.`));```

ras5_v1.mw

Edit: simpler is

`eq14 := lhs(eq13) = MM %. simplify(value(MM^(-1).rhs(eq13)));`

## convert to image...

Try

```p := plot(x^2):
im:=convert(p,Image):
ImageTools:-Embed(im);```

## dchange...

Use dchange for this. I did the equations separately, but you can put them in a list and do it all in one step.

 >
 >
 >

Desired new variables and parameters in terms of the old variables

 >

Maple needs the old variables in terms of the new variables.

 >

Try it directly

 >

We see we have to divide each side by K*R. Adding simplify simplifies the rhs, but could do this as a separate step.

 >

 >

 >

 >

 >

## sometimes possible...

`y := 2*hypergeom([1/2, -k/2], [1 - k/2], csc(omega*T)^2)`

Depending on what you know about k, you may be able to convert it to a nice form. For example, for k=3

`convert(eval(y, k = 3), StandardFunctions)`

gives

## isosurface...

implicitplot3d makes use of the ISOSURFACE structure, see ?plot,structure. I believe that this cannot be coloured with a color function, only GRID and MESH structures:

"For specification of multiple colors, the COLOR(t, A) structure, where t is RGB, HSV or HUE, is similar to that for the 2-D case, except that m by n GRID and MESH structures  are also supported. In this situation, A must be an Array with datatype=float[8] and dimensions 1..m, 1..n when t is HUE, or 1..m, 1..n, 1..3 otherwise."

Here's a workaround and how to use a colorfunction of 3 variables (simple in Maple 2016 and after)

 > restart

Desired function in implicit form

 > C := 2*x*y*z-x^2-y^2-z^2+1;

Desired colour function in 3-variable form

 > colfunc := (x, y, z) -> (x^2 + y^2 + z^2)/3:

Cannot use a colour function for implicit plot
One (ugly) workaround is to cast C in explicit form, handling the two pieces separately.
(Numeric option not required for operator plot calls as here)

 > z1,z2:=solve(C,z); f1:=unapply(z1,x,y,numeric): f2:=unapply(z2,x,y,numeric):

Find the two corresponding colorfunctions

 > colfunc1:=(x,y)->colfunc(x,y,f1(x,y)); colfunc2:=(x,y)->colfunc(x,y,f2(x,y));

 > plot3d([f1,f2],-1..1,-1..1,color=[colfunc1,colfunc2]);

For Maple 2016 and later, colorscheme can be used

 > plot3d([f1,f2],-1..1,-1..1,colorscheme=["xyzcoloring",colfunc]);

 >

## expected...

They are not to any order, but they are what I would expect. If I simplify by choosing specific examples for the initial functions, then play around with the expansion order, I find the total degrees in the residuals range from order-2 to 2*order-1. Since you have terms with a function times a first derivative of a function, you can expect up to (order plus (order-1)). For the low end, differentiation reduces the order and you have some third-order derivatives. So I would say the answer is verified up to order-3.

 > restart:
 > with(DETools):
 > PDE1 := diff(eta(t,x),t) + 1/2*diff(u(t,x),x) + 1/2*eta(t,x)*diff(u(t,x),x) - 1/48*diff(u(t,x),x\$3) + diff(eta(t,x),x)*u(t,x);

 > PDE2 := diff(u(t,x),t) + u(t,x)*diff(u(t,x),x) + diff(eta(t,x),x,t,t) + diff(eta(t,x),x) - 1/6*diff(u(t,x),x,x,t);

 > sys := rifsimp([PDE1, PDE2]);

 > id := initialdata(sys[Solved]);

Take a specific example

 > id2:=eval(eval(id),{_F1(t)=t,_F2(x)=0,_F3(x)=1,_F4(t)=1,_F5(t)=0,_F6(t)=0});

 > sols := rtaylor(sys[Solved], id2, point=[t = 0, x = 0], order = 6);

Total degrees of terms in residual, and min and max of these

 > pdetest(sols,PDE1): map(degree,[op(%)]); min(%),max(%);

 > pdetest(sols,PDE2): map(degree,[op(%)]); min(%),max(%);

 >

## userinfo...

You could use the userinfo mechanism to do this. For example,

```pkg := module() export squareme; option package; local ModuleLoad;
squareme := proc(x) x^2 end proc;
end module; ```

which leads, with infolevel[pkg] >= 2:

 > restart;
 > infolevel[pkg]:=3;

 > libname:=libname,interface(worksheetdir):
 >
 > with(pkg);

 >

## OK if c__2 depends on beta__c...

You probably know which signs of c__1, c__2 etc to choose the right solution (nicest form of solutions?).

 > restart;
 > with(DEtools):
 > ode1 := diff(c(T),T)=-2*c(T)*(1+beta__c*c(T)^2/p__c(T)^2); ode2 := diff(p__c(T),T)=2*(1+beta__c*c(T)^2/p__c(T)^2)*p__c(T); sys := [ode1,ode2]:

First solve the system with beta__c = 0.
We want the solution for non-zero beta__c to reduce to this in the limit of beta__c -> 0.

 > sys0:=eval(sys,beta__c=0); dsolve(sys0);

Now solve the system, with solving order [c(T), p__c(T)] (see ?dsolve,system), where the last one is solved first. Because we didn't use the explicit option, a partial solution is given.

To get the full solution we take any of the first solutions for p__c(T) and substitute into any of the second de's for c(T) - this is what option explicit does.

The de's for c(T) have beta__c in the denominator, which seems to be a problem, but let's proceed.

 > sol1 := dsolve(sys,[c(T),p__c(T)]);

We see that for c__2 = 0, the p__c(T) solutions reduces to the beta__c = 0 one, so we can just replace c__2 by c__2*beta (or c__2*beta^2, etc) to get solutions with the right limiting behavior.
(Perhaps c__2 = beta__c is OK? Looks simpler).

 > sol2:=[eval(dsolve(sys,[c(T),p__c(T)],explicit),c__2=c__2*beta__c)];

 > map(odetest,sol2,sys);

The denominator issue disappears for beta__c>0.

 > sol3:=simplify(sol2) assuming beta__c::positive;

Check the limiting behavior - simplify symbolic is reckless here; should consider each case for different signs of c__1, c__2 etc.

 > simplify(eval(sol3,beta__c=0),symbolic);

## twins...

I think I understand better now. You never need to remove anything, just add the latest model to the list if it doesn't match any of the others, or add it as a twin if it does. If that is correct, then this works. Since lists are inefficient, this uses a table containing Arrays for the twins and a Vector (could be an Array) for the list.

```KeyModels := Vector(): # unique models (was list l)
Models := table():     # equivalent models (was list of twins)
nkeys := 0:
for n to nmodels do
# see which key model the new model n is equivalent to
# by comparing with the existing key models
# here do it by a random number generator
r := rand(1 .. nkeys + 1)();
if r = nkeys + 1 then   # didn't match any
nkeys := nkeys + 1;
# new model - record it as key
# and start a new Array
KeyModels ,= n; # append to key models
Models[n] := Array([n]); # or perhaps Array([])
else                    # matched existing model KeyModelNum
KeyModelNum := KeyModels[r];
Models[KeyModelNum] ,= n;
end if;
end do:
```
 > restart;
 > nmodels:=100000;

 > st:=time():
KeyModels := Vector(): # unique models (was list l)
 > time()-st;

Results

 > KeyModels;

 > Models[4]; Models[28];

 >