## 60 Reputation

7 years, 154 days

## My mistake on the probability calculatio...

@vv

Yes I think it is my mistake to think that the probability of fD and fA taking the same value is  fD(x)*fA(x)dx integrated over all space. The probability isn't calculated in that way. Thank you for replying to this anyways.

## A guess at the problem...

@mapleq2013 This is a followup to the last post of the disturbing problem. I now strongly feel that it has something to do with the eigenvalue surfaces crossing each other. I originally thought that they don't cross, but since now I'm dealing with surfaces touching each other at certain locations, it's unavoidable that right before they touch they have to cross with the surfaces that are closer to each other first (which may or may not themselves touch, and may or may not at the same location). In this event it seems that I have to abandon the fast method of doing eigenvalues and sorting them. It puzzles me however that the original fsolve doesn't work either, because I think they should work by properly tracking the eigenvalues.

## Even more troublesome problems ......

I was trying to figure out why the g is not properly changing the result in one code but not in another, and I realized there is a more troublesome problem that actuallly could be the cause of this:

See in variable_not_properly_changing_results.mw , when g=0, there are 4 elements in Bm that in a particular basis (through the transformation matrices T1 and T2) are isolated that they can be extracted out to form a 2by2 matrix Br (Br for B reduced) to yield eigenvalue surfaces just for that dimension. I did it in 2by2_sample_for_g_not_properly_changing_results.mw  and it gave me the surfaces correct to my scientific intuition -- the two surfaces touch at 6 points (which is obvious in the zoomed in view). But shouldn't I be able to find this solution just by solving the original 6by6 matrix Bm for eigenvalue surfaces? The fact I couldn't find this result by setting g=0 in the original 6y6 algorithm probably means that something is wrong?

On the other hand, in the variable_properly_changing_results.mw, I originally thought this was a good code, but now I'm not so sure, because I did the same reduction to 2by2: 2by2_sample_for_g_apparently_properly_changing_results.mw and I got eigenvalue surfaces that didn't look the same as any of the 6 surfaces from solving the 6by6 problem!

These discrepancies are very disturbing because they create this contradiction between things that I believe to be correct ...

## Code modified for clarity...

@Carl Love I'm sorry about the confusion but it seems that you can Acer both didn't like it so I'll change it to this:     variable_not_properly_changing_results.mw

## Isn't C just a number times B when g=0?...

@acer  A3 is basically a number and A4 is simply zero when g=0, so C would be a number times B when g=0, am I right? I feel that this might be a stupid question, so I apologize in advance if this turns out to be due to my stupid mistake.

## g is still not changing the results...

Yes even if you put Bm in constructing C you still don't have the correct effect from g.

I'm sorry about this confusion that I put B in constructing C, I was trying to show you that even if I put B in there and picture it generates looks the same as if you put Bm in there.

## I learned so much from this, now I need ...

Wow Carl thank you so much for the in-depth explanation above! I feel that I have learned so much by reading it through, now I need some time to digest. So basically Maple is able to create an object that is not evaluated at definition, which actually feels like the function definition in C++ (I'm not a professional in coding or computer science, but I'll try to draw connection to my limited experience with C++). The things that are evaluated immediately in C++ are the inline functions I think. I have seen that people usually define C++ functions in a .h file while they call them later in a .cpp file (though I seldom do that in my scientific computing tasks because I don't feel entirely comfortable with linking files). The reason I seemed to be obsessed with loops is that in C++ looping through a variable will automatically substitute it into anything inside the loop that contains the variable, while in Maple it doesn't. So I guess in Maple when we define a procedure it's like defining it in a .h file, and when we call it later it's like in the .cpp file.

Let me digest these more.. In the meanwhile I'm trying to figure out a new problem that I encountered after using Acer's method, I'd appreciate it very much if you could have a look at it.

## A new problem arises...

Thank you for your last reply about simplifications. It's a little beyond me though because I'm currently still struggling to make things work, and Carl is writing me a deep tutorial on procedures that I try to digest ...

Here I do have a new problem that arises from using your way of solving the eigenvalues (it's super fast!), I actually think it's my lack of skill, not your way, that is causing the problem. The situation is that in the attached code, below Eq.(17), if I define C as (-1/A3).(B+A4), C should be dependent on the parameter g which is defined in the start of the code. However changing g doesn't seem to affect the surface plotted below. I'm only plotting surfaces 2 and 5 below, because I have a rough idea how they should look like. The first picture is the surfg which is surface-global, and the second is surfl which is surface-local. I know that when g=0 C is basically the same as B, and hence surfg should have those points touching each other, which should be obvious in the local view in surfl. However, even if I set g=0 or simply make C:=B, the surfaces always look the same as g not equal to zero. That's not correct because I remember earlier (before I started using your algorithms) I could do that just fine. Unfortunately I didn't save my earlier stuff because I thought your way of doing it is superior. However, I do have another code that's working on the same problem with different matrices, in there I'm using your algorithm but g is properly changing the results (for surfaces 1 and 6). I'm deeply puzzled by this inconsistency.

variable_not_properly_changing_results.mw

variable_properly_changing_results.mw

## What's the purpose of the line?...

So I guess by C := map(u-> simplify(simplify(combine(u)), size), B) you are applying a simplification to each element of B, but there are two variables in B , x and y, this map would still result in two variables? Ok I just looked at the result of the line and it has no difference to B.

But then the question becomes, what's the purpose of this line then, if it doesn't do anything apparently? Or is it the same as unapply as used by Carl to create a function out of x and y variables?

## Critical step not understood...

Thank you again for the explanations!

I just realized that there is a small but critical step in your code that I don't understand:

C := map(u-> simplify(simplify(combine(u)), size), B)

Apparently if I just do C:=B  (suppose I just use B, not  -1/A3.(B+A4) ) the loop that follows won't work. I think I don't understand either map or combine properly, and there is a size too ... Also why is there a u? What purpose does the u serve?

## Questions about how to increase calculat...

I have some questions in the following:

1. I certainly agree that if put in arrays then the computation time can be greatly reduced, but I have to admit I'm not very proficient with handling arrays in Maple, and I really want to learn what you have done there. I guess this has something to do with what I mentioned in the reply to Carl, I'm always confused with Maple's looping process and how it assigns values to arrays. e.g. in your line

for i to N do for j to N do M[i, j,  .. ] := Vector([fsolve(Dt(-Pi+2*Pi*(i-1)/(N-1), -Pi+2*Pi*(j-1)/(N-1)), w, complex)]) end do end do

you didn't call the third index of M (what does .. mean here? Isn't .. supposed to be between two numbers as a range?) or any index for the Vector, but it seems Maple automatically knows what to put in there.

2. I don't know how to properly use the matrix after it's been populated.

allReP := [seq(plots:-surfdata(Array(map(Re, M[.. ,  .. , i]), datatype = float), -Pi .. Pi, -Pi .. Pi, axes = box, color = clist[i], glossiness = .8, lightmodel = light1, style = patchcontour), i = 1 .. 6)]

This line contains a lot of things that I don't know how to use properly. So if I remove the syntax for coloring and making axes ... seq is the same as calling a loop through i=1..6? How does it know the increment is one? So you are extracting the rank 2 tensors from the rank 3 tensor M into arrays sequentially from i=1..6 and then plot with surfdata ... so what's the difference between surfdata and matrixplot? What is map? It seems to be the same as eval?

3. Yes I think sort is important here because I need the six surfaces to belong to the corresponding solutions. My eigen surfaces do not cross so a simple sort would do. I don't know why you said there is a crossing of eigenvalues since they actually repel each other with a small gap, and never crossed. However, the fact that they don't cross this time doesn't mean they won't next time, so that's why I tried to solve for it in x and y first and then plot in x and y. In this event I think maybe it's not possible to reduce computation time by using the previous method.

4. And about the floating simplify thing, I know that Maple is good for rationals and symbolics, while Matlab is for numbers. But I don't want to use Matlab for the reason I said in 3: it's hard to determine if everything in the first surface actually corresponds to the first solution if I don't first solve it in x and y. I only did the floating thing because I couldn't get the solve to work, and by changing things to floats I have observed some unreliable success in going through solve.

First thank you so much again for the detailed answers!

I have some follow-up questions:

1. "What's necessary is to compute the determinant outside the call to fsolve or else the determinant would be recomputed for every point (x,y)"

I don't quite understand what is "to compute the determinant outside the call to fsolve". My understanding is that, solve(Dt(y,x), w) would create a symbolic solution first and then substitute y and x later during the plotting process. That didn't work reliably since Maple wasn't able to find a symbolic solution everytime. I tried to use fsolve(Dt(y,x), w) myself earlier but it tries to create numerical solutions and returned me an error. I also tried to loop through all x and y first and then solve(Dt(y,x), w) which successfully calculated but then I lose track of which eigenvalue belongs to which solution (later I realized since my eigenvalues do not cross I can actually rank them from highest to lowest and got the plots, but it's very slow because Maple doesn't seem to handle massive loops as well as Matlab. So you somehow overcame all of those problems by (y,x)-> Re(fsolve(Dt(y,x), w, complex) in my understanding by letting Maple temporarily hold x and y as variables while calling fsolve and then supply x and y values later (I can feel it but would appreciate more explanations). Also wouldn't unapply be unnecessary because (y,x)-> already defines a function?

2. For the subs being necessary in solD||n:= (y,x)-> Re(fsolve(Dt(y,x), w, complex)[n])

This is my long standing confusion with using Maple to do loops. Because I have some experience progamming with C++ and in there you don't need to do subs each time (and I think in other languages you don't need to do that either). I understand Maple is not a progamming language, but what puzzles me is that how Maple doesn't substitute the n directly? Wait, sometimes it does substitute directly because I was looping through x and y and the numbers did go in solve(Dt(y,x), w). So yea I'm confused here ...

## It worked!!! Some questions about the sy...

Thank you so much Carl!! It worked!

If you don't mind I'd like to ask you some questions about the syntax you used.

1. Dt:= unapply(LinearAlgebra:-Determinant(C), (y,x)):

So the unapply here is to convert the expression in x and y into a function, is my understanding correct? If correct, why do we have to do this? I guess it has something to do with the next step but I'm not sure.

2. solD||n:= subs(_n= n, (y,x)-> Re(fsolve(Dt(y,x), w, complex)[_n]))

I don't understand why there is a _n=n, is it just to avoid using n=n? So can I use k=n alright?

Also you mentioned using a procedure, and the usually way I understand procedures is like:

function:= proc(parameter):: returnType

procedure content

end proc

But the procedure you defined is not of the above structure, so I'm a little confused.

3. The subs command, why do we need it? Can I simply write

solD||n:= (y,x)-> Re(fsolve(Dt(y,x), w, complex)[n])

Well I just tried it and it didn't work, so it seems the subs is important here, but why? And how is subs different from eval

Again thank you very much!

## @acer  Yes you were right that I go...

Yes you were right that I got that warning when I had an unassigned name in the DE without explicitly specifying it as a parameter.

Also for the 1D case, I think I did change the 2D syntax to a 1D one, but apparently Maple still considers it as 2D because I didn't convert it.

Also thank you for the explanation on the kernel and the GUI.

 1 2 Page 1 of 2