dharr

OK if c2 depends on betac...

Answered: dharr 6416

March 24 2024

2 7

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)=-2c(T)(1+beta__cc(T)^2/p__c(T)^2); ode2 := diff(p__c(T),T)=2(1+beta__cc(T)^2/p__c(T)^2)p__c(T); sys := [ode1,ode2]:

diff(c(T), T) = -2*c(T)*(1+beta__c*c(T)^2/p__c(T)^2)

diff(p__c(T), T) = 2*(1+beta__c*c(T)^2/p__c(T)^2)*p__c(T)

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

${c(T) = c__2*exp(-2*T), p__c(T) = c__1*exp(2*T)}$

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

$[{p__c(T) = -((1/2)*I)*(8*c__1*(exp(T))^8+64*c__2)^(1/4), p__c(T) = ((1/2)*I)*(8*c__1*(exp(T))^8+64*c__2)^(1/4), p__c(T) = -(1/2)*(8*c__1*(exp(T))^8+64*c__2)^(1/4), p__c(T) = (1/2)*(8*c__1*(exp(T))^8+64*c__2)^(1/4)}, {c(T) = -(1/2)*2^(1/2)*(beta__c*p__c(T)*(diff(p__c(T), T)-2*p__c(T)))^(1/2)/beta__c, c(T) = (1/2)*2^(1/2)*(beta__c*p__c(T)*(diff(p__c(T), T)-2*p__c(T)))^(1/2)/beta__c}]$

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

$[{c(T) = -(1/2)*2^(1/2)*16^(1/2)*(beta__c^2*c__2/(2*c__1*(exp(T))^8+16*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = -((1/2)*I)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}, {c(T) = (1/2)*2^(1/2)*16^(1/2)*(beta__c^2*c__2/(2*c__1*(exp(T))^8+16*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = -((1/2)*I)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}, {c(T) = -(1/2)*2^(1/2)*16^(1/2)*(beta__c^2*c__2/(2*c__1*(exp(T))^8+16*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = ((1/2)*I)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}, {c(T) = (1/2)*2^(1/2)*16^(1/2)*(beta__c^2*c__2/(2*c__1*(exp(T))^8+16*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = ((1/2)*I)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}, {c(T) = -2*2^(1/2)*(-beta__c^2*c__2/(2*c__1*(exp(T))^8+16*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = -(1/2)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}, {c(T) = (1/2)*2^(1/2)*(-32*beta__c^2*c__2/(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = -(1/2)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}, {c(T) = -2*(-2*beta__c^2*c__2/(2*c__1*(exp(T))^8+16*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = (1/2)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}, {c(T) = (1/2)*(-64*beta__c^2*c__2/(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/2))^(1/2)/beta__c, p__c(T) = (1/2)*(8*c__1*(exp(T))^8+64*c__2*beta__c)^(1/4)}]$

>	map(odetest,sol2,sys);

The denominator issue disappears for beta__c>0.

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

$[{c(T) = -2*2^(1/2)*(c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = -((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}, {c(T) = 2*2^(1/2)*(c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = -((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}, {c(T) = -2*2^(1/2)*(c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = ((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}, {c(T) = 2*2^(1/2)*(c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = ((1/2)*I)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}, {c(T) = -2*2^(1/2)*(-c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = -(1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}, {c(T) = 2*2^(1/2)*(-c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = -(1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}, {c(T) = -2*2^(1/2)*(-c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = (1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}, {c(T) = 2*2^(1/2)*(-c__2/(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/2))^(1/2), p__c(T) = (1/2)*2^(1/2)*(2*c__1*exp(8*T)+16*c__2*beta__c)^(1/4)}]$

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

$[{c(T) = -2*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = -((1/2)*I)*2^(3/4)*c__1^(1/4)*exp(2*T)}, {c(T) = 2*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = -((1/2)*I)*2^(3/4)*c__1^(1/4)*exp(2*T)}, {c(T) = -2*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = ((1/2)*I)*2^(3/4)*c__1^(1/4)*exp(2*T)}, {c(T) = 2*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = ((1/2)*I)*2^(3/4)*c__1^(1/4)*exp(2*T)}, {c(T) = -(2*I)*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = -(1/2)*2^(3/4)*c__1^(1/4)*exp(2*T)}, {c(T) = (2*I)*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = -(1/2)*2^(3/4)*c__1^(1/4)*exp(2*T)}, {c(T) = -(2*I)*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = (1/2)*2^(3/4)*c__1^(1/4)*exp(2*T)}, {c(T) = (2*I)*2^(1/4)*exp(-2*T)*c__2^(1/2)/c__1^(1/4), p__c(T) = (1/2)*2^(3/4)*c__1^(1/4)*exp(2*T)}]$

Download dsolve.mw

twins...

Answered: dharr 6416

March 22 2024

2 2

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];
        # add to existing models
        Models[KeyModelNum] ,= n;
    end if;
end do:

Download twins.mw

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

Answered: dharr 6416

March 21 2024

1 1

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

Answered: dharr 6416

March 21 2024

0 0

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

Answered: dharr 6416

March 20 2024

3 31

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

_rtable[36893490145425527140]

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

Matrix(2, 5, {(1, 1) = -1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 1, (1, 5) = 0, (2, 1) = 1, (2, 2) = -1, (2, 3) = -1, (2, 4) = 0, (2, 5) = 1})

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

[sigma__d2__new = sigma__d2/sigma__d1, sigma__v__new = gamma*sigma__d1*sigma__v, lambda__new = sigma__d1*lambda]

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

${sigma__d2__new = Psi, sigma__v__new = Gamma}$

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

[sigma__d2 = Psi*sigma__d1, sigma__v = Gamma/(gamma*sigma__d1), lambda = lambda__new/sigma__d1]

>	allvals:=[allvalues(Lambda)]: # just different index values for the RootOf

>	plot3d(allvals[1],Gamma=0..10,Psi=0..10);

>

Download nondim.mw

size...

Answered: dharr 6416

March 19 2024

1 1

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

Table of Contents...

Answered: dharr 6416

March 17 2024

4 3

You need to set up a Table of Contents to load into the help database with HelpTools:-TableOfContents:-Load(TOC,helpdbname);

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

Answered: dharr 6416

March 15 2024

1 2

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

implicitplot3d.mw

one way...

Answered: dharr 6416

March 14 2024

1 1

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

>

Add to 3:
1 1 1
Add to 4:
2 1 1
1 1 2
1 2 1
Add to 5:
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.

>

s := 0; for n to 10 do comps := Iterator:-Composition(n, parts = 3); s := s+add(1/mul(comp), `in`(comp, comps)); print(n, s) end do

>

Download compositions.mw

First one to StandardFunctions...

Answered: dharr 6416

March 13 2024

2 3

Download dchange.mw

Dirac...

Answered: dharr 6416

March 13 2024

2 2

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)=DoseDirac(t-1)-k1xSA1(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]);

>

Download Dirac.mw

Biorthogonality...

Answered: dharr 6416

March 13 2024

2 2

@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).

>

Matrix(%id = 36893490995083995364)

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)

>

revs, rvecs1 := Eigenvectors(A); evs, rperm := sort(revs, output = [sorted, permutation], key = evalf); rvecs := rvecs1[() .. (), rperm]

evs, rperm := Vector(3, {(1) = 1-sqrt(2), (2) = 1, (3) = 1+sqrt(2)}), [3, 1, 2]

Matrix(%id = 36893490995083958140)

Left (row) eigenvectors corresponding to increasing eigenvals

>

levs, lvecs1 := Eigenvectors(LinearAlgebra:-Transpose(A)); evs, lperm := sort(levs, output = [sorted, permutation], key = evalf); lvecs := lvecs1*LinearAlgebra:-Transpose([() .. (), lperm])

evs, lperm := Vector(3, {(1) = 1-sqrt(2), (2) = 1, (3) = 1+sqrt(2)}), [2, 3, 1]

Matrix(%id = 36893490995083988852)

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)

>

lvecs[1, () .. ()], rvecs[() .. (), 2]; lvecs[1, () .. ()].rvecs[() .. (), 2]; lvecs[1, () .. ()].rvecs[() .. (), 3]

Vector[row](%id = 36893490995083967172), Vector[column](%id = 36893490995083967292)

All possibilities

>

Matrix(%id = 36893490995083956340)

>

Download Biorthogonal.mw

normalization...

Answered: dharr 6416

March 12 2024

3 0

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

Vector(6, {(1) = 754.5312606638714, (2) = 512.7138273868854, (3) = 73.24600906942915, (4) = 35.82888224512065, (5) = 4.059820634692874, (6) = 0.})

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