dharr

Dr. David Harrington

8355 Reputation

22 Badges

21 years, 11 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
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.

MaplePrimes Activity


These are replies submitted by dharr

@tomleslie Yes, it is interesting how often this question is asked. Perhaps it should be in huge letters on the help pages for plot and plot3d.

@sursumCorda I had a closer look into how the algorithm works, and my analysis below suggests in the numerical case, close eigenvalues that should be detected as degenerate are not in LinearAlgebra:-MatrixFunction, but are in linalg:-matfunc. I have submitted an SCR.


For going beyond, the normal case is well-behaved and can be done just by applying the function to the eigenvalues in a unitary decomposition (Horn and Johnson, Theorem 6.2.37). So for numerical matrices with shape=symmetric or shape=hermitian (or the skew versions), that would be a separate case worth doing.

Matrix function algorithm

restart

with(LinearAlgebra)

A := Matrix([[1, 1, 3], [0, 1, 0], [0, I, 2]]); f := x/(exp(x)-1)

Matrix(%id = 36893490720140850108)

x/(exp(x)-1)

Exact answer

fA := MatrixFunction(A, f, x)

Matrix(%id = 36893490720140844196)

Matrix is non-diagonalizable, eigenvalues: 2,1,1

JordanForm(A); cp := CharacteristicPolynomial(A, x)

Matrix(%id = 36893490720140857460)

x^3-4*x^2+5*x-2

Let's try as an example first, not the function f, but a quintic polynomial

p := add((i+1)*x^i, i = 0 .. 5)

6*x^5+5*x^4+4*x^3+3*x^2+2*x+1

MatrixFunction(A, p, x); eval(p, x = A)

Matrix(3, 3, {(1, 1) = 21, (1, 2) = 70+690*I, (1, 3) = 900, (2, 1) = 0, (2, 2) = 21, (2, 3) = 0, (3, 1) = 0, (3, 2) = 300*I, (3, 3) = 321})

Matrix(%id = 36893490720205323788)

Any polynomial can be reduced to one of degree <n that gives the same result. We just find the remainder on dividing by a polynomial such that p(A)=0, so the characteristic polynomial or minimal polynomial (the same in this case) can be used

r := rem(p, cp, x); eval(r, x = A)

r := 230*x^2-390*x+181

Matrix(%id = 36893490720205308740)

mp := MinimalPolynomial(A, x)

x^3-4*x^2+5*x-2

The polynomial r is an interpolating polynomial since it is equal to p at the eigenvalues. The slopes are also equal at the doubly-degenerate eigenvalue.

plot([p, r], x = 0 .. 5, 0 .. 1000, color = [red, blue])

For a general function described by an "infinite degree" polynomial (Taylor series), we may expect it (and some slopes) to be described at the eigenvalues by a polynomial r of degree less than n. This is the interpolating polynomial.  Then f(A) = r(A)

In the present case r has to have the same values as f at the eigenvalues 1 and 2, and at the doubly degenerate eigenvalue 1 the slopes have to be the same.

f1 := eval(f, x = 1); f11 := eval(diff(f, x), x = 1); f2 := eval(f, x = 2)

1/(exp(1)-1)

1/(exp(1)-1)-exp(1)/(exp(1)-1)^2

2/(exp(2)-1)

fs := evalf([f1, f11, f2])

[.5819767070, -.3386968875, .3130352854]

Let's find a quadratic with the same values. An interpolating algorithm is used in practice, but conceptually it is simple to do it this way:

q := a*x^2+b*x+c

a*x^2+b*x+c

ans := solve({eval(q, x = 1) = f1, eval(q, x = 2) = f2, eval(diff(q, x), x = 1) = f11}, {a, b, c})

{a = -(exp(2)*exp(1)-2*(exp(1))^2-2*exp(2)+3*exp(1))/(((exp(1))^2-2*exp(1)+1)*(exp(2)-1)), b = (2*exp(2)*exp(1)-4*(exp(1))^2-5*exp(2)+6*exp(1)+1)/((exp(1)-1)^2*(exp(2)-1)), c = 2*((exp(1))^2+exp(2)-2*exp(1))/(((exp(1))^2-2*exp(1)+1)*(exp(2)-1))}

qnum := eval(q, ans); qf := evalf(%)

-(exp(2)*exp(1)-2*(exp(1))^2-2*exp(2)+3*exp(1))*x^2/(((exp(1))^2-2*exp(1)+1)*(exp(2)-1))+(2*exp(2)*exp(1)-4*(exp(1))^2-5*exp(2)+6*exp(1)+1)*x/((exp(1)-1)^2*(exp(2)-1))+2*((exp(1))^2+exp(2)-2*exp(1))/(((exp(1))^2-2*exp(1)+1)*(exp(2)-1))

0.6975546596e-1*x^2-.4782078194*x+.9904290612

eval(qf, x = 1), eval(diff(qf, x), x = 1), eval(qf, x = 2)

.5819767078, -.3386968875, .3130352862

And we get the correct exact result.

simplify(eval(qnum, x = A)); simplify(%-fA)

Matrix(3, 3, {(1, 1) = 1/(exp(1)-1), (1, 2) = ((-1+9*I)*exp(1)-(3*I)*exp(2)-1)/((exp(2)-1)*(exp(1)-1)), (1, 3) = -3/(exp(1)+1), (2, 1) = 0, (2, 2) = 1/(exp(1)-1), (2, 3) = 0, (3, 1) = 0, (3, 2) = -I/(exp(1)+1), (3, 3) = 2/(exp(2)-1)})

Matrix(%id = 36893490720140843116)

But we have to know the eigenvalues are distinct, which for numerical eigenvalues is difficult. For example, suppose we assumed the eigenvalue pair at 1 was actually two close values

eps := 0.5e-8; ff := [eval(f, x = 1.0), eval(f, x = 1+eps), eval(f, x = 2.)]

0.5e-8

[.5819767070, .5819767052, .3130352855]

Then the quadratic we find is

`~`[`=`]([eval(q, x = 1.0), eval(q, x = 1+eps), eval(q, x = 2.)], ff); ans2 := solve({%[]}, {a, b, c})

[1.00*a+1.0*b+c = .5819767070, 1.000000010*a+1.000000005*b+c = .5819767052, 4.*a+2.*b+c = .3130352855]

{a = 0.9105857850e-1, b = -.5421171570, c = 1.033035286}

qf2 := eval(q, ans2); eval(qf2, x = 1.), eval(qf2, x = 1+eps), eval(qf2, x = 2.)

0.9105857850e-1*x^2-.5421171570*x+1.033035286

.5819767075, .5819767057, .3130352860

and the matrix function is

simplify(eval(qf2, x = A))

Matrix(%id = 36893490720140862996)

Which is significantly different to the correct answer

evalf(fA)

Matrix(%id = 36893490720140849388)

In the linalg case, the input to the polyomial interpolator `matfunc/pinterp` is the list of eigevalues [1.+0.*I, 1.+0.*I, 2.+0.*I] and b, which is what I have called fs above, i.e., it is clear that the degeneracy of the eigenvalues has been detected.

On the other hand, in the LinearAlgebra case, `MatrixFunction/PolynomialInterpolation` receives the list of the eigenvalues [HFloat(0.9999999999999997)+HFloat(0.0)*I, HFloat(1.0)+HFloat(0.0)*I, HFloat(1.9999999999999998)+HFloat(0.0)*I] and the vector Vector(3, {(1) = HFloat(0.5819767068693266)+HFloat(0.0)*I, (2) = HFloat(0.5819767068693265)+HFloat(0.0)*I, (3) = HFloat(0.31303528549933135)+HFloat(0.0)*I}). The first two entries of the vector are very slightly different, suggesting that the degeneracy is not detected. (However, it could be that this detection is delegated to `MatrixFunction/PolynomialInterpolation`)

NULL

Download MatrixFunctionAlgorithm.mw

@ishan_ae I'm guessing that you didn't add two unevaluation quotes around the add: '' add '' (these are two ' not one "). This is an endless problem in Maple. One way to solve not knowing how many is just to have one in the definition, and then add a "numeric" option, which effectively tells animate not to bother unless a number is provided. The output of the fourier_f definition is now ugly.

fourier_analysis_numeric.mw

Or you can use Sum, but then you need value:

fourier_analysis_Sum.mw

 

@ishan_ae

I changed to assume n is a positive integer, then the sin(n*Pi) and cos(n*Pi) can be simplified. assume(N>0) has no effect since the N within the fourier_f is not the same as the one outside.

The delayed evaluation quotes around 'add' are because at the time the fourier_f definition is evaluated, the value of N is not known. This prevents add from generating an error because N does not yet have a value. In the animated one, I have needed to do this twice.

display is in the plots package so needs to be plots:-display (or use with(plots): at the start of the worksheet).

If you use sum instead of add, it tries to work out the general formula every time sum in encountered, which slows it down significantly.

Here it is working.

restart; k := 2; assume(n, posint)

target_f := proc (x) options operator, arrow; piecewise(-Pi < x and x < 0, 0, 0 < x and x < Pi, x, Pi < x and x < 2*Pi, 0, 2*Pi < x and x < 3*Pi, x-2*Pi) end proc

proc (x) options operator, arrow; piecewise(-Pi < x and x < 0, 0, 0 < x and x < Pi, x, Pi < x and x < 2*Pi, 0, 2*Pi < x and x < 3*Pi, x-2*Pi) end proc

target_f_plot := plot(target_f, -Pi .. 3*Pi, scaling = constrained, legend = "Target Function", color = red)

NULLa0 := simplify((int(target_f(x), x = -Pi .. Pi))/(2*Pi))

(1/4)*Pi

a_n := (int(target_f(x)*cos(n*x), x = -Pi .. Pi))/Pi

(-1+(-1)^n)/(n^2*Pi)

b_n := (int(target_f(x)*sin(n*x), x = -Pi .. Pi))/Pi

-(-1)^n/n

fourier_f := unapply(a0+(''add'')(a_n*cos(n*x)+b_n*sin(n*x), n = 1 .. N), N)

proc (N) options operator, arrow; (1/4)*Pi+('add')((-1+(-1)^n)*cos(n*x)/(n^2*Pi)-(-1)^n*sin(n*x)/n, n = 1 .. N) end proc

fourier_f_plot := plots[animate](plot, [fourier_f(N), x = -Pi .. 3*Pi, legend = "Fourier Approximation", color = blue], N = [seq(1 .. 10, 1)])

plots:-display(fourier_f_plot, target_f_plot)

``

NULL

Download fourier_analysis.mw

Not clear if this is OK - you can remove input from the worksheet using view->show/hide contents and then export to .pdf (or print to .pdf with a .pdf printer driver)

@sursumCorda Finally got a version with topological sort to work, added to the above - about twice as fast.

[Edit - I removed a redundant line of code and now can't reproduce the faster timing - it's about the same.]

Have never used these and they might not be what you want, but there is a Grading package, and Maple TA. The help page ?MapleTAIntegration has a number of useful links.

This works for me so presumably something about your installation is different. I suggest you contact technical support

compile_ex.mw

@sredh01 Your main problem is the absence of the RandomMatrix(n,n,generator=-5..5): Here's a version that generates a given number of matrices.

GivenDet2.mw

 

@ijuptilk I interpreted your question "Hi, please can someone help on how non-dimensionalize PDEs. I have tried the following, but is not working" as meaning "I know how to non-dimensionalize, and was in the process of non-dimensionalizing when Maple gave an error, so could someone please solve the Maple error". So I made the various substutions exactly as you had them. You presumably actually know the dimensions of various quantities, and I wouldn't normally attempt a non-dimensionalization without that information.

More generous MaplePrimer @sand15@mmcdara interpreted your question as requesting a complete solution to non-dimensionalize your equation. But they are slowed down in this by having to deduce things, like the f[1] and f[2] interpretation, which you have only now provided, and now we see that K[1] and K[3] appear more often than they did in the original equation. So are you OK with a solution without f[1] and f[2]? Your response is not very clear about this. You also have not commented on @mmcdara's analysis. Perhaps you now know how to proceed with the tools at hand and that is the end. But if you want more, some more information about what you know and what you want would seem appropriate. 

@ijuptilk I'm assuming the problem is in the last step. Try adding params=[K[1],mu]. I'll update the other one. Or perhaps Maple is hinting that you should take @sand15's advice :-)

@mmcdara Nice analysis. For the case where you postmultiplied the Gram Schmidt matrix by DiagonalMatrix( [1$(N-1), T*sign(d)]) you have scaled the last column by magnitude T, so they are random but with a different range from the other columns. Your next example solves that problem.

I was trying to work an integer case based on GS.DiagonalMatrix( [1$(N-1), T]).GS^+ (which solves the above problem, I think), but if T=1, then you get the identity matrix, which suddenly is not very random, even though T<>1 looks random. This made me wonder if the entries of GS itself are too correlated to be random, but that is beyond my limited statistics expertise.

@lcz Tarjan's algorithm seems to be standard, and was for directed graphs, (but works also for undirected if allowance is made for the double counting). So, yes it is strange not to have the capability for both.

It was the fastest of the algorithms I played with, but perhaps with dense graphs other algorithms might be better. Tarjans's could probably be converted to a compiled version in Maple, if the stack was implemented "by hand" using arrays. I was originally thinking algorithms based on cycle basis was the way to go, but the literature doesn't seem to point to that.

I also have another version that uses Zeon algebra, but just counts and doesn't produce the cycles themselves, so I added it to the other thread:

How-To-Determine-The-Number-Of-Cycles-Of--Graph

Infolevel gives information that is not that enlightening. 

solve.mw

I can't decide whether this is a bug, or whether we are expected to be cautious when mathematical functions like ln are involved.

First 43 44 45 46 47 48 49 Last Page 45 of 87