It isn't all clear to me what you are after. I may be misunderstanding it all the way through.
I really didn't understand what the ranges of the i and j indices of C_ik are. So I just used all eleven of them. So fix that up.
For the rest of it, well, the chances are good that I made at least some transcription error, if not worse. So look for that.
# S is a list of the first 11 nonnegative roots of -BesselJ(1,x) ?
# If so, then find these using fsolve with avoid={} option, or with plot.
S := [0., 0.3831705970, 0.7015586670, 1.017346814, 1.332369194, 1.647063084, 1.68438, 1.961585851, 2.276008438, 2.590367209, 2.904682853, 3.218967991]:
B := Vector(11,(i)->`if`(i=1,1,S[i]*coth(S[i]))):
A := Vector(11,(i)->2/(100*BesselJ(0,S[i]*10)^2)):
C := Matrix(11,11, (i,j)-> B[j]*A[i]/1 * ( evalf(Int(BesselJ(0,S[i]*x)*BesselJ(0,S[j]*x)*x^(15.0/200.0),x=0.0..0.075)) + evalf(Int(BesselJ(0,S[i]*x)*BesselJ(0,S[i]*x)*BesselJ(0,S[j]*x)*x,x=0.076..10.0)) ) );
CI := C + LinearAlgebra[IdentityMatrix](11);
LinearAlgebra[MatrixInverse](CI);
I don't see how animation comes into it.
acer
A routine which would allow efficient iteration over an rtable (Array, Matrix, Vector) as well as allow simultaneous access to the indices as each element is accessed would be very nice. It's unfortunate that it does not exist; using side-effects with rtable_scanblock is inefficient and does not scale well.
You can map or Map over an rtable, but not get access to the indices as each element is dealt with. The same goes for syntax like,
for x in A do .... od: # for rtable A
As for the NULLs in your example, select() and its friends produce the same sorts of objects that they are passed. And it seems that the size of an Array is a quality preserved by select. So, some value is needed, and NULL is better than a default of 0, since 0 might have some special meaning for a given example of selection.
acer
You can write a routine to do this, using loops. Or you could have it use rtable_scanblock, which can be made to walk over Arrays, Matrices, and Vectors in all sorts of useful ways. One neat aspect is that, while walking such a structure, rtable_scanblock has access to the indices and so can save them elsewhere as a side-effect.
For example,
FindIndex := proc(comparison,A)
local checkindex,G;
checkindex := proc(val,ind)
if comparison(val) then
G[val] := op(ind);
end if;
NULL;
end proc:
rtable_scanblock( A, [rtable_dims(A)], checkindex, [[1],A[1]]);
convert(G,list);
end proc:
foo := Array( 1..10 , i -> i ) * 2;
FindIndex( x -> evalb( x mod 3 = 0 ), foo );
You could make it more general, depending on whether you want it to handle rtables of several dimensions or to return something other than a list.
acer
Consider this sort of looping,
for i from 1 to 100 do
cat( "p", convert(5000+i,string) );
od;
So, you can produce strings, for names of files, programmatically.
acer
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[20](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[20](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
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[1]+x[2])/K) and (1-1/a) ?
L:= Matrix([[s[1]*exp((x[1]+x[2])/K)*(1-1/a),F],[s[1]*exp((x[1]+x[2])/K)*1/a,s[2]]]);
VV := Vector([x[1],x[2]]);
Z := L.VV - VV;
e1 := isolate(Z[2],x[1]);
solve(eval(Z[1],e1),{x[2]});
sol2 := op([%][2]); # pick the non-trivial solution
sol1 := op(solve(eval(Z[2],sol2),{x[1]}));
op(solve(Z[2],{x[1]}));
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[100]( eval( newZ, {K=1/2,s[1]=1,s[2]=2,a=3,F=1/2} ) );
acer
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
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
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
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[8] 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
If you use t as the dummy variable of integration alongside t[0] or t[o], then t isn't really free. This is because names like t can refer to tables in Maple, and t[0] 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
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
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
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
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[7](Xnew), results) then
results := results union {evalf[7](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