Carl Love

Carl Love

28035 Reputation

25 Badges

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

MaplePrimes Activity


These are replies submitted by Carl Love

@pik1432 

There are many ways to input a coefficient matrix using Maple's 1D input. I'd guess that there are more ways to do it than there are for inputting any other structure. In the way that I used, the rows are separated by semicolons, and the entries within each row are separated by commas. Adding line breaks and indentation is always syntactically valid, so this can be used to make the input itself look like a matrix. For example, the coefficient matrix above could be input as

C:= < 
    #xl12 xl13 xl23
      1,   1,  -1;
      1,  -1,   1;
     -1,   1,   1
>; 

As always in linear algebra, each matrix column represents the coefficients of a particular varable. 

In Matrix and Vector construction, the angle brackets ... are always used in a left and right matched pair, just like parentheses ... ), square brackets [ ... ], and curly braces ... }. Of course, these angle brackets are the same ASCII keyboard characters as the less-than and greater-than symbols.

column vector is specified in angle brackets with commas separating the entries. In the case above, I could've used

X:= <
    X1,
    X2,
    X3
>;
xl:= < 
    xl12,
    xl13,
    xl23
>;

row vector uses as the separator, for example R:= < 3 | 5 | 7 >. (The code in my Answer above didn't use any row vectors.)

A Matrix or Vector can be multiplied or divided by a scalar coefficient using the familiar and operators; so C/2 divides every entry of my above by 2.

The command LinearAlgebra:-LinearSolve(A, b) solves the matrix equation Ax = b for column vector x(It can also be used if b is a matrix.) 

The elementwise operator can be appended to another operator, such as =~, to indicate that an operation is to be applied to the elements rather than to the vectors as a whole. Scalar multiplication is automatically elementwise, but C/2 could also be written C/~2 if you want.

So, given the input above, your solution can be found as 

xl =~ LinearAlgebra:-LinearSolve(C/2, X)
          

@abdgafartunde Have you been able to turn any of this information into a procedure?

@miasene Please post a worksheet showing the "no results." Did you use the procedure Perm for the intial constructions, exactly like I had?

@Carl Love One benefit of a two-loop method is that the inner loop can accumulate the norm of the change in the X vector without needing to copy the vector. The smallness of that norm is a better convergence criterion than the residual (A.X - B). 

Thus, the body of the outer loop could be:

err:= 0; 
for i to n do 
    (X[i], t):= ((B[i] - add(A[i,j]*X[j], j= N minus {i}))/A[i,i], X[i]); 
    err:= err + abs(X[i] - t) #1-norm 
od

where N = {$1..n}.

Addendum: I delight is seeing how much a code can be compressed (still maintaining readable spacing and indentation) by using Maple 2019 or later syntax (although I think that you're using Maple 2017). Doing this, the entire procedure can be made

GaussSeidel:= proc(
    A::Matrix(square), B::Vector, X0::Vector,
    abstol::positive:= 10.^(2-Digits), maxiters::posint:= 99
)
local X:= Vector(X0, datatype= hfloat), n:= numelems(X), i, k;
    for k to maxiters do until
        abstol > 
        add(abs((X[i] - (X[i]+= (B[i] - A[i].X)/A[i,i]))), i= 1..n)
    ;
    if k > maxiters then error "did not converge" else X fi
end proc
:

The add could be parallelized to Threads:-Add if abs(...is replaced by CopySign(..., 1).

@janhardo Okay, I can read it now. Thanks. Now you pick which of those 9 curves that you want to start with and enter it into a Maple worksheet. Also enter the "other" curve that you're supposed to find intersection with. And please don't clutter up the worksheet with book pages.

@janhardo The book pages in your worksheet are too small for me to read (due to Maple not properly supporting my QuadHD screen---zooming doesn't increase the size of graphics). Please transcribe the problem that you want to work on and its instructions.

@Thomas Richard Do you know if it's possible to turn off the clickability of error messages? I hate accidentally clicking them.

@tomleslie I'm not recommending it as a workaround. I'm recommending it as the only truly correct and robust way---even if the phenomenon is deemed a bug and fixed.

I wasn't aware of this phenomenon; it is indeed shocking. However, I've argued for years that the only robust method is to always use the package prefix inside procedures. The prefix may be abbreviated, but not safely eliminated:

    uses LA= LinearAlgebra;
    LA:-MatrixInverse(M)

Also, using LA:-MatrixInverse is robust, but using LinearAlgebra[MatrixInverse] (as is recommended on nearly every help page of a package command as the "long form") is not robust.

@abdgafartunde Sorry, I was in a hurry when I wrote "I can't conceive of using two loops." Upon further reflection, it is certainly conceivable, but it is not necessary. The inner loop, which does the updating, can be replaced with a single call to LinearAlgebra:-ForwardSubstitute.

Here's an outline: Let the original problem be solving for X in A.X = B. Let X0 be an inital guess of the solution. In the procedure's initialization, split A as L + U where L is lower triangular and U is strictly upper triangular. This can be done like this:

    L:= Matrix(
        A, 'shape'= 'triangular'['lower'], 'storage'= 'triangular'['lower'],
        'datatype'= 'hfloat'
    );
    U:= Matrix(
        A, 'shape'= 'triangular'['upper'], 'storage'= 'triangular'['upper', 'strict'],
        'datatype'= 'hfloat'
    );

Now the body of the outer (and only) loop can simply be the single line

        X:= LinearAlgebra:-ForwardSubstitute(L, B - U.X)

if your convergence criterion is based on the residual A.X - B. If the criterion is based on absolute error, then you can't overwrite X immediately, and you need one more line:

        X:= copy(X1);
        X1:= LinearAlgebra:-ForwardSubstitute(L, B - U.X)

@janhardo I think that this MaplePrimes thread is getting too long. It takes a long time to load, and scrolling to the latest update is difficult. Would you please ask the above as a new Question? Indeed, I think that topically it is indeed an independent Question, worthy of its own thread. Surely we shouldn't be expected to discuss every exercise in your whole textbook in this one thread.

@janhardo In a version of this Question that you've substantially altered, you had a problem and question regarding the circle command. The reason for your error is that circle is in the plottools package, not the plots package; so, if you call it as plottools:-circle, you wouldn't have the error.

@janhardo If a curve is specified parametrically by 

x= f(t), y= g(t), -infinity < t < infinity;

then its arclength between t1 and t2 is

Int(sqrt(diff(x,t)^2 + diff(y,t)^2), t= t1..t2).

You could think of that as a path or line integral (the precise distinction between those two is inconsistently defined---it varies among authors), but I think that it's better for now if you just think of it as the 1-D integral shown above. As a VectorCalculus integral, it would be

VectorCalculus:-PathInt(1, [x,y]= 'Path'(<f,g>, t= t1..t2)).

So, when finding the arclength with this method, the "integrand" is simply 1.

@abdgafartunde Try writing some pseudocode first. Check Wikipedia.

@miasene Some response indicating whether or not my Answer was useful to you would be appreciated.

First 189 190 191 192 193 194 195 Last Page 191 of 708