If the file is already opened by Maple, then fprintf will write to its end.
restart:
fd := fopen("mres.out", WRITE):
M := Matrix([[1.,i],[-i,1]]):
with(LinearAlgebra):
for i from 1 to 10 do
res := LinearAlgebra:-Eigenvalues(M);
fprintf(fd,"%a %a \n", seq(res[j],j=1..2) );
end do:
# flush when ready, to commit the write to the file
fflush(fd);
acer
I see, thanks. So, that method may or may not work, and testing the definition (U.V1-Evals[1]*V1) is the way to find out afterwards (for each eigenvalue). It might work some or most of the time.
The first example I gave seems to be in the problematic space.
restart:
with(LinearAlgebra):
U := Matrix([[a^2+b^2-1,1,0],[1-a^2-b^2,1,0],[1,0,1]]):
Evals := simplify(Eigenvalues(U), {a^2 + b^2 = 1}):
U3:= U - Evals[3]*IdentityMatrix(3):
V3:= LinearSolve(U3[1..2,1..3],Vector(2)):
V3:= eval(V3,_t[3]=1); # or whatever the parameter is.
[0]
[ ]
V3 := [0]
[ ]
[1]
map(simplify,U.V3-Evals[3]*V3,{a^2+b^2=1}); # not OK
[0]
[ ]
[0]
[ ]
[1]
Using the side relations during pivot selection seems to work. I'd advise trying it that way.
Normalizer:=x -> simplify(simplify(x,{a^2+b^2=1})):
LinearSolve(U3,Vector(3));
[-_t0[3]]
[ ]
[ 0 ]
[ ]
[_t0[3] ]
V3:= eval(%,_t0[3]=1); # or whatever the parameter is.
[-1]
[ ]
V3 := [ 0]
[ ]
[ 1]
map(simplify,U.V3-Evals[3]*V3,{a^2+b^2=1}); # OK
[0]
[ ]
[0]
[ ]
[0]
acer
I see, thanks. So, that method may or may not work, and testing the definition (U.V1-Evals[1]*V1) is the way to find out afterwards (for each eigenvalue). It might work some or most of the time.
The first example I gave seems to be in the problematic space.
restart:
with(LinearAlgebra):
U := Matrix([[a^2+b^2-1,1,0],[1-a^2-b^2,1,0],[1,0,1]]):
Evals := simplify(Eigenvalues(U), {a^2 + b^2 = 1}):
U3:= U - Evals[3]*IdentityMatrix(3):
V3:= LinearSolve(U3[1..2,1..3],Vector(2)):
V3:= eval(V3,_t[3]=1); # or whatever the parameter is.
[0]
[ ]
V3 := [0]
[ ]
[1]
map(simplify,U.V3-Evals[3]*V3,{a^2+b^2=1}); # not OK
[0]
[ ]
[0]
[ ]
[1]
Using the side relations during pivot selection seems to work. I'd advise trying it that way.
Normalizer:=x -> simplify(simplify(x,{a^2+b^2=1})):
LinearSolve(U3,Vector(3));
[-_t0[3]]
[ ]
[ 0 ]
[ ]
[_t0[3] ]
V3:= eval(%,_t0[3]=1); # or whatever the parameter is.
[-1]
[ ]
V3 := [ 0]
[ ]
[ 1]
map(simplify,U.V3-Evals[3]*V3,{a^2+b^2=1}); # OK
[0]
[ ]
[0]
[ ]
[0]
acer
There's no provision in this for avoiding using hidden zeros as pivots, during the linear system solving. I still thnik that, in general, that's a serious issue for this class of problem.
Did you mean to have,
V1:= LinearSolve(U1[1..3,1..4],Vector(3));
but perhaps got caught by the < issue?
This seems like just a bit more "by-hand" work, than calling NullSpace.
Also, how can you know that it's "OK" to drop the last row of U1 (as opposed to the 2nd, etc)? What if the 1st and 2nd rows are linear combinations of each other and the 4th is independent of the top 3? Why not just do,
LinearSolve(U1,Vector(4));
which is in effect what NullSpace does?
acer
There's no provision in this for avoiding using hidden zeros as pivots, during the linear system solving. I still thnik that, in general, that's a serious issue for this class of problem.
Did you mean to have,
V1:= LinearSolve(U1[1..3,1..4],Vector(3));
but perhaps got caught by the < issue?
This seems like just a bit more "by-hand" work, than calling NullSpace.
Also, how can you know that it's "OK" to drop the last row of U1 (as opposed to the 2nd, etc)? What if the 1st and 2nd rows are linear combinations of each other and the 4th is independent of the top 3? Why not just do,
LinearSolve(U1,Vector(4));
which is in effect what NullSpace does?
acer
Interesting. Using a simplifier (with side relations) for Normalizer on your Matrix U resulted in a spurious all-zero Matrix of eigenvectors when using the Eigenvectors command. So maybe Eigenvectors is not clever enough about how it goes to work in such a situation.
But I still think that you might be rightfully concerned with Eigenvectors() using hidden zeros as pivots during the kernel computation.
Having computed the eigenvalues as evU, you might then compute the eigenvectors associated with the i'th eigenvalues the good old fashioned way. That is to say, using NullSpace().
You may also find that the Normalizer setting gets the eigenvalues in their simplified form too.
with(LinearAlgebra):
U := Matrix(1/8*[[0, 0, -8*I, 0], [0, 8*I, 0, 8*I], [I*b, I*(a-1), (2*I)*(a-1+b)+8*I, (2*I)*(a-1-b)], [I*(a-1), -I*b, (2*I)*(a-1-b), -(2*I)*(a-1+b)]]):
use_assump := x -> simplify(simplify(x, {a^2+b^2 = 1})):
Normalizer := use_assump:
evU := Eigenvalues(U);
for i from 1 to 4 do
NullSpace(U-evU[i]*IdentityMatrix(4))
end do;
Normalizer := normal:
acer
Interesting. Using a simplifier (with side relations) for Normalizer on your Matrix U resulted in a spurious all-zero Matrix of eigenvectors when using the Eigenvectors command. So maybe Eigenvectors is not clever enough about how it goes to work in such a situation.
But I still think that you might be rightfully concerned with Eigenvectors() using hidden zeros as pivots during the kernel computation.
Having computed the eigenvalues as evU, you might then compute the eigenvectors associated with the i'th eigenvalues the good old fashioned way. That is to say, using NullSpace().
You may also find that the Normalizer setting gets the eigenvalues in their simplified form too.
with(LinearAlgebra):
U := Matrix(1/8*[[0, 0, -8*I, 0], [0, 8*I, 0, 8*I], [I*b, I*(a-1), (2*I)*(a-1+b)+8*I, (2*I)*(a-1-b)], [I*(a-1), -I*b, (2*I)*(a-1-b), -(2*I)*(a-1+b)]]):
use_assump := x -> simplify(simplify(x, {a^2+b^2 = 1})):
Normalizer := use_assump:
evU := Eigenvalues(U);
for i from 1 to 4 do
NullSpace(U-evU[i]*IdentityMatrix(4))
end do;
Normalizer := normal:
acer
Since nobody has yet stated it in just so many words, I will add that `=` does not by itself bring about assignment in Maple.
So, it is not true that, "The symbol = means assignment." Without qualification an assertion like that is completely general. But since it isn't true in Maple in particular, then it also isn't always true.
Assignments in Maple are done using a `:=` statement, or using the `assign` operator.
As for what `=` can in fact mean in Maple, please see other replies to this subthread whose points I won't repeat here.
acer
Since nobody has yet stated it in just so many words, I will add that `=` does not by itself bring about assignment in Maple.
So, it is not true that, "The symbol = means assignment." Without qualification an assertion like that is completely general. But since it isn't true in Maple in particular, then it also isn't always true.
Assignments in Maple are done using a `:=` statement, or using the `assign` operator.
As for what `=` can in fact mean in Maple, please see other replies to this subthread whose points I won't repeat here.
acer
The procedure Rsquared below could allow you to compute R^2 for data appearing in the first and second columns of a Matrix and for a given fitting formula in variable t.
Rsquared := proc( formula, T::symbol, m::Matrix )
local N, ybar, SS_err, SS_tot, SS_reg;
Digits:=trunc(evalhf(Digits));
N := op([1,1],m);
ybar := add( m[i,2], i=1..N )/N;
SS_tot := add( ( m[i,2] - ybar )^2, i=1..N );
SS_reg := add( ( eval(formula,T=m[i,1]) - ybar )^2, i=1..N );
SS_err := add( ( m[i,2] - eval(formula,T=m[i,1]) )^2, i=1..N );
return 1-SS_err/SS_tot;
end proc:
M := Matrix(6,2, [[1,2],[2,3],[3,4.8],[5,10.2],[6,30.9]] );
form := 0.888 + 0.606*exp(0.649*t);
Rsquared( form, t, M );
You could fix it up to do proper argument checking and validation. Comments could be added to the code. Or you could just examine it and compare to the formulas on that web-page, and try and figure out a bit more of how the Maple programming language works.
acer
The procedure Rsquared below could allow you to compute R^2 for data appearing in the first and second columns of a Matrix and for a given fitting formula in variable t.
Rsquared := proc( formula, T::symbol, m::Matrix )
local N, ybar, SS_err, SS_tot, SS_reg;
Digits:=trunc(evalhf(Digits));
N := op([1,1],m);
ybar := add( m[i,2], i=1..N )/N;
SS_tot := add( ( m[i,2] - ybar )^2, i=1..N );
SS_reg := add( ( eval(formula,T=m[i,1]) - ybar )^2, i=1..N );
SS_err := add( ( m[i,2] - eval(formula,T=m[i,1]) )^2, i=1..N );
return 1-SS_err/SS_tot;
end proc:
M := Matrix(6,2, [[1,2],[2,3],[3,4.8],[5,10.2],[6,30.9]] );
form := 0.888 + 0.606*exp(0.649*t);
Rsquared( form, t, M );
You could fix it up to do proper argument checking and validation. Comments could be added to the code. Or you could just examine it and compare to the formulas on that web-page, and try and figure out a bit more of how the Maple programming language works.
acer
So, you may be interested in the so-called "coefficient of determination", or R^2.
See
this link for details and explanations.
What I showed before is the residual sum of squares, which that page describes as its SS_err term. That pages also gives some formulae for R^2. Perhaps you would like to use the formula which computes R^2 via regression sum of squares.
If your X and Y data lie in the first and second columns of Matrix M, then in your code you can access the entries as M[i,1] instead of X[i] and M[i,2] instead of Y[i]. That saves your having to form the X and Y Vectors explicitly.
acer
So, you may be interested in the so-called "coefficient of determination", or R^2.
See
this link for details and explanations.
What I showed before is the residual sum of squares, which that page describes as its SS_err term. That pages also gives some formulae for R^2. Perhaps you would like to use the formula which computes R^2 via regression sum of squares.
If your X and Y data lie in the first and second columns of Matrix M, then in your code you can access the entries as M[i,1] instead of X[i] and M[i,2] instead of Y[i]. That saves your having to form the X and Y Vectors explicitly.
acer
The uploaded file has,
p*n+s = `mod`(0, a)
Did you perhaps intend either of these forms below?
> restart:
> s := 3: p := 1: a := 4: t := 0:
> for n to 10 do
> if `mod`(p*n+s,a) = 0
> then t := t+1
> else p*n+s
> end if;
> end do;
> t;
3
Or,
> restart:
> s := 3: p := 1: a := 4: t := 0:
> for n to 10 do
> if p*n+s mod a = 0
> then t := t+1
> else p*n+s
> end if;
> end do;
> t;
3
acer
Exactly, Scott. It's not currently possible to get the symbol to display correctly in both places; Subject line and Recent comments sidebar. It's nice to know that Will's looking at some issues with this site.
As for the motivation, I was actually trying to represent a key press in a not so unusual way, eg. -. But it's also pretty natural to imagine that the less-than symbol might be wanted, at some point, in the title of a posting on a discussion site for mathematical software.
acer
