## 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.

## no eps for 3D plots...

Maple does not properly export .eps for 3D plots (works fine with vector graphics for 2D plots). You can export .pdf from the plot context menu and it will work directly in overleaf with the graphicx package. If you want to do it programmatically you can use the jpeg or png driver, but of course these are bitmap formats (as is the .pdf)..

## most have equivalents...

Most things in linalg have their equivalent in LinearAlgebra, so I wouldn't say it is split in two. Note that the ArrayTools package is restricted to rectangular storage array objects (at least in 2023). For these, ArrayTools:-IsZero is fast.

For sparse (or rectangular) matrices, M-> andmap(`=`,M,0) is a simple alternative.

[Edit: see @sursumCorda's response - andmap is only good for moderately sized rectangular matrices]

## non-dimensionalize...

The non-dim Lambda is a function of 2 non-dim variables, so you will need 3D plots this time.

[Edit in blue - the choice of variables to eliminate was not working; was always the first 2]

 > restart;
 > local gamma:with(LinearAlgebra):

Want to "non-dimensionalize" this expression.

 > p3 := RootOf((64*gamma^2*sigma__d1^2*sigma__d2^6*sigma__v^2 + 64*gamma^2*sigma__d2^8*sigma__v^2 + 512*sigma__d2^6)*_Z^10 - 2*gamma^6*sigma__d1^2*sigma__v^6 - 8*gamma^5*sigma__d1^2*sigma__v^5*_Z + (-3*gamma^6*sigma__d1^4*sigma__v^6 - 3*gamma^6*sigma__d1^2*sigma__d2^2*sigma__v^6 + 24*gamma^4*sigma__d2^2*sigma__v^4)*_Z^2 + (-12*gamma^5*sigma__d1^4*sigma__v^5 - 52*gamma^5*sigma__d1^2*sigma__d2^2*sigma__v^5 - 48*gamma^5*sigma__d2^4*sigma__v^5 + 32*gamma^3*sigma__d1^2*sigma__v^3 + 160*gamma^3*sigma__d2^2*sigma__v^3)*_Z^3 + (12*gamma^6*sigma__d1^4*sigma__d2^2*sigma__v^6 + 36*gamma^6*sigma__d1^2*sigma__d2^4*sigma__v^6 + 24*gamma^6*sigma__d2^6*sigma__v^6 - 12*gamma^4*sigma__d1^4*sigma__v^4 - 216*gamma^4*sigma__d1^2*sigma__d2^2*sigma__v^4 - 412*gamma^4*sigma__d2^4*sigma__v^4 + 32*gamma^2*sigma__d1^2*sigma__v^2 + 384*gamma^2*sigma__d2^2*sigma__v^2)*_Z^4 + (32*gamma^5*sigma__d1^4*sigma__d2^2*sigma__v^5 + 224*gamma^5*sigma__d1^2*sigma__d2^4*sigma__v^5 + 248*gamma^5*sigma__d2^6*sigma__v^5 - 336*gamma^3*sigma__d1^2*sigma__d2^2*sigma__v^3 - 1392*gamma^3*sigma__d2^4*sigma__v^3 + 384*gamma*sigma__d2^2*sigma__v)*_Z^5 + (4*gamma^6*sigma__d1^6*sigma__d2^2*sigma__v^6 + 16*gamma^6*sigma__d1^4*sigma__d2^4*sigma__v^6 + 16*gamma^6*sigma__d1^2*sigma__d2^6*sigma__v^6 + 4*gamma^6*sigma__d2^8*sigma__v^6 + 16*gamma^4*sigma__d1^4*sigma__d2^2*sigma__v^4 + 512*gamma^4*sigma__d1^2*sigma__d2^4*sigma__v^4 + 1072*gamma^4*sigma__d2^6*sigma__v^4 - 176*gamma^2*sigma__d1^2*sigma__d2^2*sigma__v^2 - 2288*gamma^2*sigma__d2^4*sigma__v^2 + 128*sigma__d2^2)*_Z^6 + (48*gamma^5*sigma__d1^4*sigma__d2^4*sigma__v^5 + 96*gamma^5*sigma__d1^2*sigma__d2^6*sigma__v^5 + 32*gamma^5*sigma__d2^8*sigma__v^5 + 512*gamma^3*sigma__d1^2*sigma__d2^4*sigma__v^3 + 2464*gamma^3*sigma__d2^6*sigma__v^3 - 1792*gamma*sigma__d2^4*sigma__v)*_Z^7 + (32*gamma^4*sigma__d1^4*sigma__d2^4*sigma__v^4 + 208*gamma^4*sigma__d1^2*sigma__d2^6*sigma__v^4 + 96*gamma^4*sigma__d2^8*sigma__v^4 + 192*gamma^2*sigma__d1^2*sigma__d2^4*sigma__v^2 + 3136*gamma^2*sigma__d2^6*sigma__v^2 - 512*sigma__d2^4)*_Z^8 + (192*gamma^3*sigma__d1^2*sigma__d2^6*sigma__v^3 + 128*gamma^3*sigma__d2^8*sigma__v^3 + 2048*gamma*sigma__d2^6*sigma__v)*_Z^9):

First write the polynomial with variable lambda, and normalize so one of the terms is 1 (arbitrarily take the first one)

 > p:=expand(subs(_Z=lambda,op(p3))): pnorm:=expand(%/op(1,%));

So all the other terms must now be dimensionless. We use the automated non-dimensionalization procedure of E. Hubert, G, Labahn, Found Comput Math 13 (2013) 479–516, doi:10.1007/s10208-013-9165-9. See also https://youtu.be/Nl2FBAbU1pE

 > terms:=[op(2..-1,pnorm)]: allvars:=[indets(p3,name)[],lambda]; #lambda as a true variable is last. A simpler result may occur if the ones to eliminate are first nvars:=nops(allvars); nterms:=nops(terms)

Matrix K contains the exponents of the variables in each term. It has 5 rows but rank 3, so we can eliminate 5 - 3 = 2 variables.

By default these are the first 2 in allvars, so rearrange allvars if you want a different result.

 > K := Matrix(nvars, nterms, (i, j) -> diff(terms[j], allvars[i])*allvars[i]/terms[j]); r1 := Rank(K);

Hubert's matrix A has full row rank, equal to the number of variables to eliminate

 > H, U := HermiteForm(K, output = ['H', 'U'], method = 'integer'): Determinant(U): A := U[-nvars + r1 .. -1, .. ]; A . K: # should be zero rA := Rank(A);

Calculate matrix W. Procedure ColumnHermiteForm is in startup code

 > HA, V := ColumnHermiteForm(A): W := V^(-1):
 > Vi := V[.., 1 .. rA]: Vn := V[.., -nvars + rA .. -1]: Wu := W[1 .. rA, ..]: Wd := W[-nvars + rA .. -1, ..]:
 > nondims := table(): nondimvars := table(): for i to nvars - rA do     nondimvars[i] := cat(allvars[i + rA], __new);     nondims[i] := nondimvars[i] = mul(`^`~(allvars, V[ .. , i + rA])); end do: nondimvars := convert(nondimvars, list); nondims := convert(nondims, list);

 > rewrites := table(): for i to nvars do     rewrites[i] := allvars[i] = mul(`^`~(nondimvars, Wd[ .. , i])); end do: rewrites := convert(rewrites, list);

So the non-dimensionalized p is

 > newsys := subs(rewrites, p);

The naming here is not optimal. sigma__v__new is similar to Gamma previously but with sigma__d1 instead of sigma__d, so Gamma or maybe Gamma_1; lambda__new should be Lamba, and the ratio sigma__d2/sigma__d1 could be some other capital Greek letter, say Psi
So the new Lambda is

 > newnames:={sigma__d2__new = Psi, sigma__v__new=Gamma}; Lambda := RootOf(eval(newsys,newnames),lambda__new);

Find how to write the old variables in terms of the new ones, and directly use eval as a check. We only need to eval for the variables that we are not eliminating; the others will automatically disappear

 > map(lhs,remove(x->rhs(x)=1,rewrites)): OldToNew := eval(solve(nondims,%)[],newnames); sigma__d1^4*simplify(eval(p, OldToNew)): # may need to scale to get rid of last of the old variables RootOf(%,lambda__new): %-Lambda; # should be zero

 > allvals:=[allvalues(Lambda)]: # just different index values for the RootOf
 > plot3d(allvals[1],Gamma=0..10,Psi=0..10);

 >

## size...

The default appears to be 500x500 pixels. If you add the plot option size=[1000,1000] then the exported file will be 1000x1000 pixels.

In the TOC, you set a priority for each entry, so make your overview the highest priority

See here on Mapleprimes for an example.

## implicitplot3d...

Since you only know M via the equation and not explicitly, implicitplot3d is appropriate.

implicitplot3d.mw

## one way...

Here's a way that only generates valid possibilities for i, j, k. I'm assuming each of i, j, k is in the range 1..n.

 >
 >

Compositions of k with 3 parts (sum to k). Need ones that add to k=3, 4, .. n

 >

Check possibilities

 >

1 1 1
2 1 1
1 1 2
1 2 1
3 1 1
2 1 2
2 2 1
1 1 3
1 2 2
1 3 1

Calculate the required sum

 >

So work this out for the first few n.
If we are going to do all of them we can avoid duplicate iterators and incrementally accumulate sums.

 >

 >

## First one to StandardFunctions...

 > restart;
 > q:=hypergeom([1, 1, 3/2], [5/2, 5/2, 3, 3], -(x/2)**2);

Note sure about the assumptions - choose one root.
Often conversions work when the argument is simple.

 > itr:=u=-x^2/4; tr:=x=solve(u=-x^2/4,x)[1];

 > qu:=PDEtools:-dchange(tr,q);

 > convert(qu, StandardFunctions); q2:=PDEtools:-dchange(itr,%,simplify);

 >

## Dirac...

If I understand what you want, Dirac can do this for you. Note that aside frim the infinitesimal time interval aspect @C_R mentioned, the ifelse immediately finds that t is not 1 and so Inf =0 after the ifelse statement is executed.

 > restart;
 > k12:=0.0452821; k21:=0.0641682; k1x:=0.00426118; Dose:=484; #Inj:=ifelse(t=1,Dose,0); ode_sys:=diff(SA2(t),t)=SA1(t)*k12-SA2(t)*k21,          diff(SA1(t),t)=Dose*Dirac(t-1)-k1x*SA1(t); ics:=SA2(0)=0,SA1(0)=0;

 > Solution1:=dsolve({ode_sys,ics},{SA1(t),SA2(t)},type=numeric);

 > plots:-odeplot(Solution1,[[t,SA1(t),color="red",thickness=1],[t,SA2(t),color="blue",thickness=1]],t=0..50,labels=["Time t","Salicylic acid"],labeldirections=[horizontal,vertical]);

 >

## Biorthogonality...

@mehdi jafari You asked about the relationship between left and right eigenvectors, and getting a diagonal matrix from them. They are biorthogonal, meaning that left eigenvectors are orthogonal to right eigenvectors corresponding to different eigenvalues. So the product of the matrices of left eigenvectors (as row vectors) and right eigenvectors (as column vectors) is diagonal.

[Edit - 2023 version uploaded] Rerun the worksheet in your version of Maple; the rtable output from Maple 2024 doesn't display correctly in earlier versions (or here on Mapleprimes, SCR submitted).

 >
 >
 >

Right (column) eigenvectors corresponding to increasing eigenvalues (any order for the eigenvalues will do, so long as it is consistent for left and right eigenvectors)

 >

Left (row) eigenvectors corresponding to increasing eigenvals

 >

Biorthogonality means that left eigenvectors are orthogonal to right eigenvectors corresponding to different eigenvalues.

e.g, left eigenvector for 1 is orthogonal to right eigenvalues for 1-srqt(2) and 1+sqrt(2)

 >

All possibilities

 >

 >

## normalization...

The eigenvectors from Eigenvectors are arbitrarily scaled. In this case they happen to be normalized but differ in signs. Once you account for this (and a forgotten square root), there is agreement. Notice you are also relying on the orders of the eigenvalues and singular values being consistent.

WHY DOES THE TRANSPOSE OF Vt DETERMINED BY THE FUNCTION SingularValues NOT AGREE WITH THE evectors CALCULATED BY THE FUNCTION Eigenvectors?

 > restart: with(LinearAlgebra):
 > X := Matrix([[8.79,9.93,9.83,5.45,3.16],            [6.11,6.91,5.04,-0.27,7.98],            [-9.15,-7.93,4.86,4.85,3.01],            [9.57,1.64,8.83,0.74,5.80],            [-3.49,4.02,9.80,10.00,4.27],            [9.84,0.15,-8.99,-6.02,-5.31]],            datatype=float[8],order=Fortran_order);
 (1)
 > U,S,Vt:= SingularValues(X, output = ['U', 'S', 'Vt'],thin=true)
 (2)
 > SDM:= DiagonalMatrix(S[1..5],5,5)
 (3)

THIS EQUALS TO ORIGINAL X MATRIX

 > U.SDM.Vt
 (4)
 > X -~ U.SDM.Vt
 (5)

EIGENVALUES

 > S*~S
 (6)

EIGENVECTORS

 > Transpose(Vt)
 (7)

COMPARE LINE (6) AND THE evalues OF  LINE  (8), IN AGREEMENT.

COMPARE LINE (7) AND THE evectors OF LINE (8), NOT IN AGREEMENT - differ in sign

 > evalues, evectors:= Eigenvectors(Transpose(X).X)
 (8)

Find signs of first rows of evectors and Transpose(Vt)

 > signsev:=convert(map(signum,simplify(fnormal(evectors[1,..]),zero)),list); signsV :=convert(map(signum,simplify(fnormal(Vt[..,1]),zero)),list); diffsigns:=zip(`*`,signsev,signsV);
 (9)
 > Transpose(Vt)-evectors.DiagonalMatrix(diffsigns);
 (10)
 > sdm:= DiagonalMatrix(sqrt~(evalues[1..5]),5,5)
 (11)

THIS SHOULD EQUAL TO THE ORIGINAL X MATRIX

 > U.sdm.Transpose(evectors.DiagonalMatrix(diffsigns));
 (12)
 > X - U.sdm.Transpose(evectors.DiagonalMatrix(diffsigns));;
 (13)
 >

## some examples...

Not sure about documentation or the level of explanation you want, but if examples help, here are some:

 >

infix -> prefix

 >

 >

 >

 >

 >

 >

postfix -> prefix

 >

 >

 >

 >

 >

more...

 >

 >

 >

 >

 >

 >

 >

## spacecurve...

A workaround is to use spacecurve instead of plot3d. Unfortunately, the curves need to be done separately.

Perpendicular_3D_lines.mw

## 3D PDF from Maple 3D plots...

Interesting; I didn't know about such a thing. I looked into it and played around a bit, and various of the exported plot formats can be used to make a u3d file, which can then be imported into Acrobat to make a 3D PDF or used in pdflatex to make a 3D PDF.

Here are the steps I tried

1. Make a 3d plot of a red cube, and export it to a file

```p:=plots:-display(plottools:-cuboid([0,0,0],[1,1,1],color=red));
Export("redcube.off", p, base=currentdir);```

I tried .off, .ply, .stl via the Export command and .x3d geometry and COLLADA (.dae) via the plot context menu. (These formats are all known to MeshLab.)

2. Convert to .u3d

I used MeshLab from www.meshlab.net (Windows 64 bit but other OSs are available). File -> Import mesh -> choose file. .ply gave an error, .stl asked about removing duplicate vertices, which I accepted. Only .off showed as a red cube; the others had no color.

3a. Open Acrobat Pro and create a blank page .pdf (using file->create). Find the rich media tool and open it. Select add 3d and drag a rectangle where you want the object to be. Choose the .u3d file. You will have to overide a media safety warning, then you can click on the image to activate it. Then click and you can move around. There was no color, and some faces of the cube were transparent. (A sphere from complexplot3d also lost color but was faithfully reproduced.)

3b. The .u3d file can instead be incorporated into .pdf via the LaTeX package media9 and pdfLaTeX. I used overleaf (www.overleaf.com), which uses TeX Live. The .zip file here contains the .tex file, the .u3d file and the generated .pdf file.

3D_PDF.zip

Comments - only the 3D objects are exported, e.g. from POLYGONS PLOT3D structures. These can be colored face by face, but by default the cuboid command sets them all at once to red. I tried coloring different faces differently, but Maple was not exporting the colors to the .obj file (and maybe others).

## no value...

The default value is shown when there is no argument passed corresponding to that parameter. (Using NULL achieves the same thing because the NULL is removed from the sequence.)

 >
 >

 >

 >

 >

 >

 >