## 17428 Reputation

14 years, 141 days

On google groups. (comp.soft-sys.math.maple and sci.math.symbolic)

On stackoverflow.

On math.stackexchange.com.

## could you be more specific?...

It's not clear whether you expect to create elements which have 20 decimal places of information, or whether you just want the Matrix to be able to hold such numbers. Matrix(2,2,(i,j)->evalf(1/(i+j+1)),datatype=sfloat); If the datatype of the Matrix is 'sfloat', which stands for software float, then the entries can have any number of digits. The precision of what can be assigned into the entries, after creation of such a Matrix, would depend on the environment variable Digits. See the help-page, ?UseHardwareFloats . Other ways to get a datatype=sfloat Matrix are, restart: UseHardwareFloats:=false: Matrix(2,2,(i,j)->evalf(1/(i+j+1)),datatype=float); MatrixOptions(%,'datatype'); restart: Digits := 20: Matrix(2,2,(i,j)->evalf(1/(i+j+1)),datatype=float); MatrixOptions(%,'datatype'); acer

## Oh, L is L(x)...

Sorry, I don't know why, but I missed that L is L(x) a function of vector x. Let's have another shot at it, then. There seems to be a missing multiplication sign, in your definition of L, in the [1,1] entry. Should there be a * between exp((x+x)/K) and (1-1/a) ? L:= Matrix([[s*exp((x+x)/K)*(1-1/a),F],[s*exp((x+x)/K)*1/a,s]]); VV := Vector([x,x]); Z := L.VV - VV; e1 := isolate(Z,x); solve(eval(Z,e1),{x}); sol2 := op([%]); # pick the non-trivial solution sol1 := op(solve(eval(Z,sol2),{x})); op(solve(Z,{x})); newZ := eval(Z,{sol1,sol2}); # Now, how to simplify that to zero? # Are there some assumptions on the other unknowns, to give # simplify() a helping hand? # This is crude. evalf( eval( newZ, {K=1/2,s=1,s=2,a=3,F=1/2} ) ); acer

## Nullspace, and Rank...

If L-IdentityMatrix(n) has full rank and if its determinant is not zero, then yes, only the trivial solution of the zero-Vector should get returned. I suspect that, in such a case, there is no other solution to be found, by any method. But if L-IdentityMatrix(n) is not of full rank and has a determinant of zero, then NullSpace should give a basis for the solution space, or LinearSolve (or `solve` of your fully formed explicit system) should give parametrized solutions. In other words, I don't think that the method is what's important here. What should be key is whether there exist non-trivial solutions (or not). acer

## ImportMatrix...

Alternatively, without having to muck about with the mouse, or assistants, MM := ImportMatrix("4399_p0000000.xls",datatype=anything); Of course, you may need to make the string point directly to the .xls file with its full path, or to change currentdir() before making the call. acer

## rearrange...

Does this follow? x = L . x 0 = L . x - Identity . x (L-Identity) . x = 0 So, in Maple commands, with(LinearAlgebra): X := LinearSolve( L-IdentityMatrix(n), ZeroMatrix(n) ); where n is ColumnDimension(L). Or you could get a basis for the nullspace in which x lives, with, NullSpace( L-IdentityMatrix(n) ); acer

## ImportData assistant...

This should work in Maple 10, as well as 11. Open the ImportData assistant, either through the Standard GUI's menu bar (Tools -> Assistants -> Import Data), or using the explicit command ImportData() . After doing that, I selected your .xls file. I changed the choice for the datatype from float to 'anything', so that it could hold the first row of names. Upon return from this assistant, I had a 2001x8 Matrix, which I assigned by issuing, M := %; acer

## because t is the dummy variable of integ...

If you use t as the dummy variable of integration alongside t or t[o], then t isn't really free. This is because names like t can refer to tables in Maple, and t references the 0 index entry of t. Hence mixing t and t[...] is a muddle to be avoided, with semantic baggage that confuses Maple. acer

## storage suitable for pointplot3d...

On reflection, it may not be at all obvious to the new user how to store the data and the final iterates, so that they can be easily plotted. restart: N := proc(f) local Xnew, X0, incX, A, B, k, df, count, Lreal, Lcplx, R; R := Array(1..101*101,1..3); df := diff(f, x); # These two outermost loops are used to create the X0 initial points. count := 0; for A from 0 by 0.02 to 2 do for B from -1 by 0.02 to 1 do X0 := evalf(A + B*I); # Now do ten Newton iterations, using X0 Xnew := X0; for k to 10 do incX := evalf(eval(f/df, x = Xnew)); Xnew := Xnew - incX; end do; count := count + 1; R[count,1],R[count,2],R[count,3]:=evalf(A),evalf(B),Xnew; end do; end do; # Return the results in lists suitable for pointplot3d. Lreal := [seq([R[k,1],R[k,2],Re(R[k,3])],k=1..101*101)]; Lcplx := [seq([R[k,1],R[k,2],Im(R[k,3])],k=1..101*101)]; return Lreal,Lcplx; end proc: Lreal,Lcplx := N(sin(x)): plots[pointplot3d](Lreal,axes=boxed); plots[pointplot3d](Lcplx,axes=boxed); acer

## what is the given function...

So, you had to use Newton's method to get from X0 to X1, X2, etc? If so, then what was the function? You must know, since you already solved question a). Someone else here also asked for this. No doubt it was laid out earlier, perhaps in Part I or Part II. Question b) seems to be that you repeat the task like in a), for each X0=A+B*I, but a double-loop doing that was already illustrated here. You might try studying it some more. acer

## save all the final iterates...

So, you want to save all the final (10th) iterates, for each A and B? If that's true, then simply create an Array, outside the double loops. (Is it 101x101 in size?) Then after each k-loop finishes, save the final Xnew to the A,B coordinate of the Array. Then return that Array at the end of the procedure. This Array would then contain the 10th and final Newton iterates, for all the X0=A+B*I initial points. You could then do some sort of complex plot. Perhaps the idea was to show you how the various subregions of the A+B*I complex domain gets mapped to final iterates. ps. For anyone who was shocked that the derivative of f is computed each time through the triple loop, sorry! Of course, it would be so much more efficient to get diff(f,x) just once, outside all the loops. For shame! acer

## newton's method, with multiple initial p...

I'm guessing that the idea is to start Newton's method using each of the points X0=A+B*I in the complex plane. So, the chances are greater to find a point that converges to a root. And perhaps several roots may be found. N := proc(f) local Xnew, X0, incX, A, B, k, results; results := {}; # These two outermost loops are used to create the X0 initial points. for A from 0 by 0.02 to 2 do for B from -1 by 0.02 to 1 do X0 := evalf(A + B*I); # Now do ten Newton iterations, using X0 for k to 10 do incX := evalf(eval(f, x = X0))/evalf(eval(diff(f, x), x = X0)); Xnew := X0 - incX; X0 := Xnew; end do; # If it's "good enough", and not found already, # then put the final iterate in the set of results. if abs(incX) < 0.1000000000*10^(-6) and not member(evalf(Xnew), results) then results := results union {evalf(Xnew)}; # Print which initial point converged to which new result. printf("initial point %Ze converged to %Ze\n", evalf(A + B*I), Xnew); end if; end do; end do; # Return the results return results; end proc: # try it out N(sin(x)); Look, this may not be what you were assigned to do. It was unclear. If it's close, then study it, then when you understand it you should be able to alter it appropriately. acer

## not a condition number issue...

The difference between the exact and floating-point eigenvector results, for say your 5x5 example, is not a condition number issue. The reason for the differences you see in the exact vs floating-point results is simply that, for any given eigenvalue, the eigenvectors live in a linear subspace. Multiples of single eigenvectors, or linear combinations of eigenvectors in that relevant subspace (eigenspace) also happen to be valid eigenvectors for that given eigenvalue. Take your 5x5 exact example. The exact eigenvalues are {3/4,0,23/20,1,1}. The eigenvalues {3/4,0,23/20} each have a single eigenvector associated with them, in the corresponding column of result Matrix e. You can scale each of those, and the result will still be an eigenvector of the same eigenvalue. The proof is easy, and falls out of the linear quality of the definition, M * x = lambda * x Similarly for the double eigenvalue {1,1}. It has two linearly independent eigenvectors associated with it (not always true, but true in this example). Any linear combination of such basis eigenvectors will produce another eigenvector in the same eigenspace associated with that same eigenvalue (ie, 1). This too falls right out of the linear quality of the definition above. I don't think that numerical conditioning has much, if anything, to do with your 5x5 example when computed in floating-point. I use Linux, and the particular scaling or linear combinations done on Windows might differ. But here is a sequence of elementary column operations with which I was able to transform the eigenvectors of evalf(L) into the floating-point equivalent of the results of the exact Matrix L. Lf := evalf(L): vf, ef := Eigenvectors(Lf): Digits:=15: QQ:=ColumnOperation(simplify(ef),1,1/0.4472135955): ColumnOperation(QQ,3,1/0.5920027785,inplace): ColumnOperation(QQ,4,-1/(.7071067812),inplace): ColumnOperation(QQ,2,-1/0.8164965809,inplace): ColumnOperation(QQ,[5,2],0.03887187786,inplace): ColumnOperation(QQ,5,(0.5/0.706304985706392),inplace): ColumnOperation(QQ,[2,5],-1,inplace): ColumnOperation(QQ,5,1/0.5,inplace): ColumnOperation(QQ,[5,2],1,inplace):Digits:=10:map(fnormal,QQ); Note that I only had to scale the 1st, 3rd, and 4th columns, relating to the eigenvalues of multiplicity one, to get an accord. To handle the eigenvalue 1, which has multiplicity two, I took linear combinations of just those the 2nd and 5th columns. You may also enjoy playing with crude forward error estimates, using the definition, ie, Norm( Lf . ef - ef . DiagonalMatrix(vf) ); You may also enjoy playing with the EigenConditionNumbers routine. These might help, if you doubt the floating-point calculations for your 50x50 case. acer

## WinXP for Dell PCs...

An article about Dell's offering WinXP for new purchases is here. acer

## data is not floating-point...

After I run your worksheet, the data in L is exact numeric, not floating-point. If floating-point results are acceptable to you, then just do something like Eigenvectors(evalf(L)) . If you really need exact results, then consider that the characteristic polynomial of L is of degree 50. This polynomial, whose roots are the eigenvalues of L, does not factor quickly, if at all. So it may take a long time trying to do so, if tried for exact L. Suppose that the characteristic polynomial doesn't factor; would you then be satisfied with a nullspace (basis of eigenspace) which was rife with unindexed implicit RootOfs? What could one do usefully with such a beast (other than to evalf it, in which case evalf(L) could have been used instead). acer

## what type of data?...

Is your data symbolic or exact rational, or numeric and floating-point, or...? It can make a difference with regard to what sort of suggestions might be useful. You might also consider posting or uploading some example. acer
﻿