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

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

 >

Notes: In principle, you could skip the KeyModels Vector and just use for i,_ in eval(Models) do to iterate through the keys (but not in numerical order). I've anticipated you might want the key model in with the twins, but if not replace Array([n]) with Array([]).

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

 >

 >

 >

 >

 >

 >

 >