dharr

Dr. David Harrington

9127 Reputation

22 Badges

21 years, 241 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

@Harry Garst So if I have some Matrices, I might want to manipulate them with their entries. In your Zehfuss example, I can see why you might want to use that rule to evaluate Determinant(Kronecker(A,B) since using that rule gives a much simpler expression than the direct evaluation. But for that purpose, you could write a procedure DetKron:=(A,B)->Determinant(A)^etc.

I can see a different application where you want to simplify some expression with abstract matrices A,B that are just the symbols A,B but have the properties of matrices (in particular non-cummutativity) - such as implementing these Knonecker properties.

You seem to want to make your life difficult by wanting to have A,B real matrices and then not use that fact - requiring uneval and uneval quotes. So everything has to be done inside a call to ZehfussReplace or similar. Or you can use inert operators and then use value on your result. But it is awkward and not the way I tend to use Maple.

For the abstract case, you can do some things with define or pattern matching. But I think the physics package is a better option. QuantumOperators act like matrices, and I suspect that Knonecker products are also there somewhere since in an abstract sense they are just tensor products. So these rules may already exist. I don't use the package much so I can't be more specific.

In general (for symbols or real matrices), if I want to know if two things are equal, I simplify their difference and see if it is zero (or the zero matrix). And in general I manipulate things in abstract form and when I get to the end of those manipulations I put in the numbers with eval. That is a bit harder with matrices, but try

A . B;
eval(%, {A = Matrix(2, 2, symbol = a), B = Matrix(2, 2, symbol = b)});

or for an expanded polynomial

2*x^3 + x^2 + 2;
eval(%, x = Matrix(2, 2, [1, 0, 3, 5]));

@acer That's much better. I struggled with getting the two square matrices into the same type specification. For the inerts, I was thinking that the OP might be wanting an inert output, e.g.,

ZehfussReplace := proc(expr) uses LinearAlgebra;
  subsindets(expr,
    specfunc(And(specfunc('Matrix'('square'),{%KroneckerProduct,KroneckerProduct}),
                 2 &under nops),{%Determinant,Determinant}),
    proc(q) local A:=op([1,1],q), B:=op([1,2],q);
      %Determinant(A)^(upperbound(B)[1])*%Determinant(B)^(upperbound(A)[1])
    end proc);
end proc:

and then value could be used on the output.

@Harry Garst Yes you can have such a type, as I illustrate in the answer I posted.

@Alfred_F Here's a relatively simple way to get a reasonable estimate. Interesting that many of the parameters are zero.

NLPSolveMatrixForm2.mw

@Harry Garst You need an extra eval in the evalb. I guess verify does that internally.

There is no point doing expr:= something; and then ZehfussReplace(expr) since expr is already evaluated. So you have to put the expression you want directly in the ZehfussReplace call - see the example in the attached.

Maple has the usual scoping rules for nested procedures - the inner ones can see the variables in the more outer ones. But the question you asked has more to do with subsindets, which scans through an expression and applies the inner procedure whenever it comes across "Determinant"

Dharr_reply2.mw

@Alfred_F You want Maximize not maximize to use the routine in the Optimization package. Your -50 in the constraint should be outside the integral, which solves the issue that @sand15 was having.

The integrals make the problem much more difficult. You have three ways to use NLPSolve (which is what Maximize uses if the problem is nonlinear) - enter the objective and constraints as expressions (your method), enter as procedures, or use the Matrix Form. I could only get it to work in the expression form if the integrals were solvable symbolically. The procedure form incorrectly counts the dummy integration variable as a problem variable, which leaves the more complicated MatrixForm. I use NLPSolve a lot, but not with constraints and I didn't get the constraintjacobian feature to work. In the end I couldn't get your very complicated expression to work. One could implement the third argument "needc" feature in the procedures, but it wasn't immediately obvious to me how to proceed.

It is worth formulating the problem with all parameters positive if you can and then use assume = nonnegative.

I think you also want to formulate the problem so your function wn comes down to zero at xend. Otherwise the integrals will capture a part of the function going below the axis. Perhaps piecewise could solve this.

I think your nested procedures confuse NLPSolve's attempts to process the objective and constraints, so I put everything as just simple expressions.

Here's some notes on how to use the various modes. Suggest you work up in complexity. 7 variables in a very nonlinear problem that isn't nondimensionalized is a very difficult problem. I'd guess the presence of the exp is a problem.

NLPSolveMatrixForm.mw

@Harry Garst By Maple's full evaluation rule, if A and B are already matrices, then after

expr := Determinant(KroneckerProduct(A,B));

expr is a number with the Determinant of the Kronecker product already done, i.e., it is too late to process it.

So the processing routine has to use ZehfussReplace := proc(expr::uneval)

With that modification your ZehfussReplace routine works at least for the simple case of providing

ZehfussReplace(Determinant(KroneckerProduct(A, B)))

It could be made more robust, e.g. checking for square matrices etc, only two matrices, allowing for symbolic A,B etc. Not sure what your precise requirements are.

If I convert the first line to 1D math you see it is

dim := 1;
A := Matrix(dim, dim, symbol = a);
evalb(factor(Determinant(`⨂`(A, A))) - abs(A)^(2*RowDimension(A)) = 0);

so I think the |A| in the 2D math is being interpreted as abs and not Determinant the first time. Somehow if you leave it as 2D, it interprets as abs the first time and then if you reexecute the same line it "has now figured out" that A is a Matrix and you get the determinant. Very mysterious. The solution is to write out Determinant(A).

for dim = 4, the Kronecker product is 16x16 and the determinant has 16! = 20922789888000 terms, and then you are trying ro factor the result so it isn't unexpected that it takes too long.

I'm not clear about what you mean by "detect whether there is a determinant of a Kronecker product present in an expression."

@Paras31 The Maple values are not what I am getting; I am using default epsilon. 

wval = 0.0008418042;

dwval = -0.062794051; (same as MMa)

theata := 0.2715071460;

FokasTests.mw

Perhaps I am not understanding the fourier plot, which looks like the MMa plot. It has a point at x = 0.0786577523618090, theta = 0.241280506060024 on the 1 min curve. This seems closer to the Maple result theta = 0.2715071460 than the MMa result 0.355. Are you sure the MMa number is correct?

The estimated error suggests that 0.0008418042 should have most of its digits correct. It is possible because of inaccurate evaluation of the integrand, that the estimated error is not correct, which might happen because of the very complicated expression. But I get the same result after normal(V) or collect(V,exp) which are quite different. I also tried getting rid of all preprocessing by making a procedure that does the Re directly in a numerical procedure; presumably what MMa does, and got the same result.

I also tried splitting the two integrals differently (different s ranges, and other corresponding changes in the lambda procedures), but also got exactly the same result for wval in each case.

Notice that the two integrals have similar magnitudes but opposite signs, which leads to larger errors in their difference wval. Then in the ratio dwval/wval you are dividing a large number by an error-prone much smaller value, which is why the error in theta is amplified.

Yes, I miss those indicators, too. My solution is to move to the edit tab and watch the "mode" when I am using mixed text/non-executable math. Then when I am entering text I see "text"; I use F5 to change to non-executable math and I see the change to "non-executable"; two times F5 goes back to text and I see "text" again. Edit: here's what I mean:

Edit: I just noticed the small icon on the drop-down box on the right of the top "line" of the screen also changes as you "toggle" F5.

I hadn't noticed differences in the cursor.

@Ronan It's a while since I used them; I can't remember :-(

@Paras31 Those seem to have good agreement, enough that you wouldn't see the differences on a plot. Your plots are different orientations and scales so hard for me to see where you think there is a difference, aside from the points at x=8 cm, which being on the boundary are exceptional. 

Are there places where the MMa and Maple results are very different but not on the boundary? Perhaps x=7.9 and t= 1min?

@Paras31 The approx_w(0.01,0.1) agrees to 3 decimal places with MMa, so is not a good sample point for further tests. Just let me know a point where there is poor agreement between MMa and Maple, and I'll take a closer look.

@Paras31 Your routine gets 0.3752737697 instead of 0.375077. Your s range is different so this is not an exact comparison, but in any case are you saying you are sure MMa is more accurate?

@Paras31 It is hard to say since I haven't followed in detail the two methods. I assume they are both approximations, so do you have some reason to think one is closer to the real solution? Have you benchmarked for a system for which you know the real solution? Can you improve the approximation, e.g., why is one integral from -5..5? is that adjustable?

First of all, I think you need to do an exact comparison (overlaid plots) or accurate comparisons at selcted points. Is there good accuracy for some x/t regimes?

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