Dr. David Harrington

4791 Reputation

21 Badges

19 years, 33 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 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.

MaplePrimes Activity

These are replies submitted by dharr

@philroe For run to run consistency you can use sort on the list. It's not always clear what order you get, but for the same list elements your should get the same order each time.

@Gharhoud That's correct. (As in @vv's first response.)

This has been improved in Maple 2023, where ExtremePoints(f(x)); returns [-1,0,2]


So this is just multiplication of each term by exp(k*2*Pi*I*(degree of that term)). So the procedure (checked against your outputs) is below. So here k=0 would give sheet1 unchanged, k=1 would give sheet2, k=2 would give sheet3. 

  local var,deg;
  # find indeterminate
  if nops(var)<>1 then
      error "unable to determine variable"
  end if;
  # deg finds degree of a term
  if generator::`+` then
    generator*exp(k*2*Pi*I*deg(generator)) # if only one term
  end if;
end proc:


You could just pass the variable through as a second argument, which will be more efficient if you are doing it many times.

@2cUniverse The aspect ratio does seem to be slightly different from the one exported with the context menu.

Edit: I added scaling=constrained to the display and then both the context menu export and the exportplot command gave the same aspect ratio.

@22117147 Happy to help. I added the new conjecture at the end. Keeping solutions in symbolic form removes the doubt that may exist with numerical solutions.


@NIMA112 That's why it is a reply and not an answer. As I said "I don't know precisely what you want, since you do not show any calculation that didn't work." Your title mentions a numerical error but your worksheet does not show any explicit error. What is your question? What exactly do you think the "singularity problem" is? I am guessing that the physical problem might have a singularity and therefore no amount of mathematics will make it go away.

I don't know precisely what you want, since you do not show any calculation that didn't work. But visually there does seem to be a singularity in the physical problem at (-1,1.5) where the contour lines converge. So I guess you can go closer by using a higher setting of Digits. Note that you have mispelled Digits so it is not actually at 30.

@Nicole Sharp I don't use the Maple add-on in Excel so don't have an answer for all the inconsistencies you see. So just two comments that might be helpful

Since Excel works in hardware precision (about 15 digits), if Maple passes it a result with 32 digit precision Excel must truncate that in further calculations, so how the extra ones are displayed or if they are just removed is a moot point.

9*10^9 is exact and not a floating point number with finite precision, so displays as "9000000000". evalf converts it to a floating point number. But 9.*10^9 has an explicit decimal point and so is a floating point number. In Maple, numbers created with the "e" or "E" exponent notation are always floating point, even if there is no explicit decimal point.

@SHIVAS As in my last reply, pdsolve can't find an exact solution, so I don't know how to find Theta_exact. Do you have some reason to believe there is an analytical solution? If so, you can perhaps use hints or other methods to help pdsolve it, but you would need to know or guess something about the solution.


For your new equation,

pdsolve(OdeSys, Theta(X, tau))

gives no solution so it is unlikely that one exists.


I put almost all the code in the startup code region, and in the click, changed etc code regions I just call a procedure defined in the startup region. This means you aren't debugging multiple bits of code in multple locations, and you can have extensive code without any problem.

I loaded the packages in the startup region and it seemed to work consistently.


@mmcdara Interesting approach. In the small amount of reading I did, I recall something about the fact that values need not be near the corresponding eigenvalue. For contour plots, I'm not sure the results are invalid, just that you might need very fine countours to catch the unusual locations?

@rzarouf These sorts of full evaluation problems for making plots are common and are best solved by using the procedure form of the plot call, though you also found a bug in Norm, for which I will submit an SCR.



A := LinearAlgebra:-DiagonalMatrix([-1., I + 1., -I]);

Matrix(3, 3, {(1, 1) = -1., (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1.+1.*I, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -I})

After replacing the regular inverse with the matrix inverse, we still have a problem

implicitplot((0.1)^(-1) < Norm(MatrixInverse((x + y*I)*IdentityMatrix(3) - A), 2), x = -1 .. 1, y = -1 .. 1);

Error, (in factor/remember) floating point coefficients are incompatible with field extension; use 'real' or 'complex' instead

Track it down

B:=MatrixInverse((x + y*I)*IdentityMatrix(3) - A);

Matrix(3, 3, {(1, 1) = 1/(x+I*y+1.), (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1/(x+I*y+(-1.-1.*I)), (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1/(x+I*y+I)})

Next is a bug in Norm


Error, (in factor/remember) floating point coefficients are incompatible with field extension; use 'real' or 'complex' instead

But if we have values for x and y it works



Aside from the bug, we don't want to do the inverse symbolically and then put numbers in (except for very small matrices perhaps). We solve this problem, really a full evaluation issue, by making a procedure.

f:=(x,y)->Norm(MatrixInverse((x + y*I)*IdentityMatrix(3) - A), 2)-(0.1)^(-1);

proc (x, y) options operator, arrow; LinearAlgebra:-Norm(LinearAlgebra:-MatrixInverse((x+I*y)*LinearAlgebra:-IdentityMatrix(3)-A), 2)-0.1e2 end proc

implicitplot(f, -2 .. 2, -2 .. 2,axes=boxed);

Doing it by singular values is more efficient; in general evaluation of matrix inverse should be avoided. The 2-norm of the resolvent (z*I-A)^(-1) is the reciprocal of the minimum singular value of z*I-A. I didn't track the following error down, but again we should avoid full evaluation by making a procedure.

implicitplot((0.1)^(-1) < 1/min(SingularValues((x + y*I)*IdentityMatrix(3) - A)), x = -1 .. 1, y = -1 .. 1);

Error, (in LinearAlgebra:-Eigenvalues) expecting either Matrices of rationals, rational functions, radical functions, algebraic numbers, or algebraic functions, or Matrices of complex(numeric) values

fs:= (x,y)->1/min(SingularValues((x + y*I)*IdentityMatrix(3) - A))-(0.1)^(-1);

proc (x, y) options operator, arrow; 1/min(LinearAlgebra:-SingularValues((x+I*y)*LinearAlgebra:-IdentityMatrix(3)-A))-0.1e2 end proc

implicitplot(fs, -2 .. 2, -2 .. 2,axes=boxed);

We could avoid some reciprocals by doing

fs2:= (x,y)->min(SingularValues((x + y*I)*IdentityMatrix(3) - A))-(0.1);

proc (x, y) options operator, arrow; min(LinearAlgebra:-SingularValues((x+I*y)*LinearAlgebra:-IdentityMatrix(3)-A))-.1 end proc

implicitplot(fs2, -2 .. 2, -2 .. 2,axes=boxed);


Download pseudospectrum2.mw

@Carl Love Thanks. Sometimes I hope the OPs provide at least some information to help. (If it's not in Horn and Johnson it seems it must be a bit exotic.) Turns out it was easier than I thought to calculate.

1 2 3 4 5 6 7 Last Page 1 of 45