Carl Love

Carl Love

28050 Reputation

25 Badges

12 years, 349 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

Here's a crude procedure for it:


restart:

PowerMethod:= proc(
     A::Matrix(square, complexcons),
     {tolerance:= 10^(1-Digits)},
     {maxiterations:= 99},
     $
)
uses LA= LinearAlgebra;
local
     n:= LA:-RowDimension(A),
     b:= Vector(n, fill= 1.),
     T,
     e0:= 0,
     e1
;
     to maxiterations do
          T:= A.b;
          e1:= LA:-Norm(T);
          if abs(e1-e0) < tolerance then  return e1  end if;
          e0:= e1;
          userinfo(2, PowerMethod, e0);
          b:= T/e0
     end do;
     FAIL
end proc:           

 

A:= LinearAlgebra:-RandomMatrix(3, shape= symmetric, datatype= float[8]);

A := Matrix(3, 3, {(1, 1) = 67., (1, 2) = -31., (1, 3) = 92., (2, 1) = -31., (2, 2) = 44., (2, 3) = 29., (3, 1) = 92., (3, 2) = 29., (3, 3) = 99.})

PowerMethod(A);

176.421771487843387

LinearAlgebra:-Eigenvalues(A);

Vector(3, {(1) = -33.5577318719757, (2) = 67.1359603841324, (3) = 176.421771487843})

 


Download PowerMethod.mw

You can do it with LinearAlgebra:-QRDecomposition or LinearAlgebra:-LinearSolve(..., method= QR).

There is some bug here; I don't know specifically where. The bug is manifested regardless of whether there are assumptions.

Here's a workaround: convert to a list and then back to a Vector.

k:= Vector[column]([W, (1/3)*W, (1/3)*W, (1/3)*W]);
Vector(sort(convert(k, list), (a,b)->is(a>b)));

The Power method won't converge for this matrix because there are two eigenvalues with the maximal absolute value.

You should stop using the deprecated linalg package. Everything you need is in LinearAlgebra, for example Eigenvalues and Eigenvectors.

You just have some minor syntax problems:

  1. You need an extra pair of square brackets to enclose the data.
  2. You used two independent variables, x and v, in the curve.

CurveFitting:-LeastSquares(
     [
          [1.0, 2.33], [2.0, 0.626e-1], [3.0, -2.16],
          [4.0, -2.45], [5.0, -.357], [6.0, 2.21],
          [7.0, 2.75], [8.0, .636], [9.0, -2.45]
     ],
     x,
     curve = a+b*cos(x)+c*sin(x)+d*cos(2*x)+e*sin(2*x)
);

Note that this solution does not explicitly use the LinearAlgebra package. It can be modified to use LinearAlgebra:-LeastSquares if that's what you want.

You should stop using the deprecated linalg package in all your code.


restart:

Pts:= [[2,0], [6,1], [8,0]]:

n:= nops(Pts):

V:= [a, b, c]:

f:= unapply(add(V[k]*x^(n-k), k= 1..n), x);

proc (x) options operator, arrow; a*x^2+b*x+c end proc

Eqns:= (p-> f(p[1]) = p[2]) ~ (Pts);

[4*a+2*b+c = 0, 36*a+6*b+c = 1, 64*a+8*b+c = 0]

A:= LinearAlgebra:-GenerateMatrix(Eqns, V, augmented);

A := Matrix(3, 4, {(1, 1) = 4, (1, 2) = 2, (1, 3) = 1, (1, 4) = 0, (2, 1) = 36, (2, 2) = 6, (2, 3) = 1, (2, 4) = 1, (3, 1) = 64, (3, 2) = 8, (3, 3) = 1, (3, 4) = 0})

R:= LinearAlgebra:-ReducedRowEchelonForm(A);

R := Matrix(3, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = -1/8, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 5/4, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1, (3, 4) = -2})

Sol:= convert(V =~ R[..,-1], list);

[a = -1/8, b = 5/4, c = -2]

eval(f(x), Sol);

-(1/8)*x^2+(5/4)*x-2

 


Download Interpolate.mw

I used exact arithmetic in your worksheet, and I got all zeros at the end. To use exact arithmetic, I did three things:

  1. I changed all numbers of the form 210e9 to 210*10^9
  2. I changed all numbers of the form 0.5, 0.25, etc., to 1/2, 1/4, etc.
  3. I removed all evalf.

Luckily, there was no significant increase in the computation time. With these changes, the setting of Digits is insignificant.

Here is the modified worksheet: Exact_arithmetic.mws

 

a,b,c:= 'RootOf(2012*x^3+2013*x+2014, index= k)' $ k= 1..3:
evala((a+b)^3+(b+c)^3+(c+a)^3);
                              3021
                              ----
                              1006


a:=["0101101","0000101","0001001"];

["0101101", "0000101", "0001001"]

A:= Matrix(length(a[1]), nops(a), (i,j)-> parse(a[j][i]));

A := Matrix(7, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (2, 1) = 1, (2, 2) = 0, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (4, 1) = 1, (4, 2) = 0, (4, 3) = 1, (5, 1) = 1, (5, 2) = 1, (5, 3) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (7, 1) = 1, (7, 2) = 1, (7, 3) = 1})

a||(1..3):= seq(A[..,j], j= 1..3);

a1, a2, a3 := Vector(7, {(1) = 0, (2) = 1, (3) = 0, (4) = 1, (5) = 1, (6) = 0, (7) = 1}), Vector(7, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 1, (6) = 0, (7) = 1}), Vector(7, {(1) = 0, (2) = 0, (3) = 0, (4) = 1, (5) = 0, (6) = 0, (7) = 1})

a1+a2+a3;

Vector(7, {(1) = 0, (2) = 1, (3) = 0, (4) = 2, (5) = 2, (6) = 0, (7) = 3})

 


Download String_to_Matrix.mw

The package is intended to be used in conjunction with Matlab. The name of every procedure in the package is the same as a Matlab command. However, you can use the package on its own, with no connection to Matlab. The main benefit to that (that I can see) is that many of the calling sequences are significantly shorter. For example, MTM:-eig is much less to type than LinearAlgebra:-Eigenvectors.

L:= ["1000","1","1110"]:

map(x-> cat("0" $ 7-length(x), x), L);

Lookup "Shoelace formula" in Wikipedia. This will give the area of any polygon as long as the vertices are given in order (clockwise or counterclockwise).

Area:= proc(P::list([realcons,realcons]))
local k, n:= nops(P);
     abs(add(P[k,1]*P[k+1,2] - P[k+1,1]*P[k][2], k= 1..n-1) + P[n,1]*P[1,2] - P[1,1]*P[n,2])/2
end proc:

Area([[0,0], [10,0], [10,10]]);
                               50

If you want to work within the geometry package, you could do this:

AreaOfTriangle:= proc(P::list([realcons,realcons]))
uses G= geometry;
local A,B,C,ABC;
     G:-area(G:-triangle(ABC, [G:-point(A, P[1][]), G:-point(B, P[2][]), G:-point(C, P[3][])]))
end proc:

AreaOfTriangle([[0,0], [10,0], [10,10]]);
                               50

There is no direct command that I know of, but the computation is trivial. The following procedure returns a sequence of the [center, radius] pairs. Isn't that good enough for you?

Gershgorin:= proc(A::Matrix(square))
local n:= LinearAlgebra:-RowDimension(A), i, j;
     seq([A[i,i], add(abs(A[i,j]), j= 1..n) - abs(A[i,i])], i= 1..n)
end proc:
Gershgorin(LinearAlgebra:-RandomMatrix(4));
          [48, 186], [-12, 149], [-12, 52], [-30, 110]

The LinearAlgebra:-Eigenvectors command returns the matrix P and the diagonal of D.

restart:
macro(LA= LinearAlgebra):
A:= LA:-RandomMatrix(4, datatype= float[8]):
E,P:= LA:-Eigenvectors(A):
Need to make D local because it is reserved for differentiation.
local D:= LA:-DiagonalMatrix(E);
LA:-Norm(P.D.P^(-1) - A);

Essentially 0, but there is some round-off error.

Use the ?Grid package rather than the Threads package. Then there are no worries about "thread safety" because it is not a shared-memory environment. In particular, use ?Grid,Map . Example:

Integrands:= [seq(sin(k*x), k= 1..2^9)]:
Grid:-Map(evalf@Int, Integrands, x= 0..Pi);

Using this, I got 100% processor utilization on my 8-CPU machine.

First 340 341 342 343 344 345 346 Last Page 342 of 395