## dharr

Dr. David Harrington

## 6416 Reputation

20 years, 23 days
University of Victoria
Professor or university staff

## Social Networks and Content at Maplesoft.com

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

## rho and rho__v...

I changed the rho's in the solution to rho__v's (as in the Eqs) and it worked.
solution_check.mw

## Solve for three variables....

I agree with Carl. Perhaps the following helps to illustrate the point.

 >
 >

 >

 >

We want to find when Q = N for any u. Since doubling the a's can be compensated for by halving the b's we can fix the scale by arbitrarily choosing a value for b[0], say 1. Or we can just write the solution for the other variables in terms of b[0]. Since we have 3 equations in 6 unknowns, we can extend this logic to another two variables aside from b[0], i.e., as @Carl Love noted we need to choose three variables to solve for in terms of the other three. For example, the following gives a[1], a[2] and b[2] in terms of the other 3. (Solve/identity is just a shortcut for the 3 equations equating the coefficients.

 >

Suppose we give a 4th variable to solve for - now we see the solution contains b[0] = b[0], meaning we can choose b[0] arbitrarily, the others are as above.

 >

In this case, Maple chose to make a[1] arbitrary; we get the same as if we had chosen a[2], b[1] and b[2] to solve for.

 >

 >

## 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);