Carl Love

Carl Love

28025 Reputation

25 Badges

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

MaplePrimes Activity


These are replies submitted by Carl Love

@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.

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