Carl Love

Carl Love

28025 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

fsolve needs a bit of help here. I'm not sure why; it may be a precision issue. I didn't investigate because the workaround was obvious to me. Define A = exp(-2*a). This gives the following system of three equations, which fsolve solves immediately:

fsolve({c*A^4+18 = 36, c*A+18 = 55, A= exp(-2*a)}, {A, a, c});
             {A = 0.7864846673, a = 0.1200910258, c = 47.04478236}

Do

Sinus:= evalf@Units:-Standard:-sin:

Then compare:

Sinus(30*Unit(degree));
Sinus(30); #30 radians

Maple is a "symbolic computation system", so those string manipulations are totally unnecessary; just use subs instead.

Your equation can be easily solved with all variables symbolic and using a logarithm to an arbitrary base, B:

restart:
equation1:= filler(a*x+b) - filler(c*x+d) = rs:
equationAct:= subs(filler= log[B], equation1):
simplify(solve(equationAct, x));
             (b - B^rs*d)/(c*B^rs - a)

The condition for both log arguments to be positive is (b*c - a*d)/(c*B^rs - a) > 0.

If you ever again get the urge to use string manipulations to do math, ask here for an easier way.

Your procedure has highly symbolic components such as isolate; it could never beat pure integer arithmetic.

Here's my version of repeated squaring. It's like VV's except that I unrolled his recursion into a do loop. It works for any associative operator (such as square matrix multiplication or modular multiplication), with the default being ordinary `*` multiplication.

RepSq:= proc(X, N::nonnegint, `*`:= :-`*`, Identity:= 1)
local x:= X, R:= Identity, n:= N;   
     do
          if n::odd then R:= `*`(x,R); n:= n-1 fi;
          if n=0 then return R else n:= n/2 fi;
          x:= `*`(x,x)
     od
end proc
: 

 

Some pointers on algorithm timing, referring to mistakes that you made in yours:

1. Don't display output and make a timing at the same time. Suppress the output and display it later if it's needed. For example, note that the number needs to be converted to base 10 to be displayed, and you don't want to be timing that.

2. You need much, much larger examples and many, many more of them to get meaningful timings. It's inconceivable that your puny numbers took 0.04 seconds to exponentiate; all you're measuring is overhead. Try a list of a thousand numbers with exponents in the hundreds of thousands. (Make absolutely sure that you suppress the output in this case!! When the GUI gets its output buffer overloaded, there's nothing that you can do except either wait hours for the output to display or kill the GUI (which unfortunately kills all your open worksheets).)

3. CPU time is a more meaningful comparison than real time for raw algorithms, such as this. Real time is more meaningful for comparing parallelized code or code that does disk I/O operations.

No Maple assumptions are needed for this:

ysol:= piecewise(x < 1/2, 2, x >= 1/2, 0); #1/2, NOT 0.5
​​​​​​diff(ysol, x);

In Maple, you should never use a decimal number such as 0.5 as an abbreviation for an exact fraction such as 1/2

Your code had some minor syntax errors, as corrected by Tom. Your algorithm (as opposed to your code) has two major flaws (not errors): It uses matrix inversion, and it unnecessarily recomputes the symbolic Jacobian on each iteration. Here is my rewrite that avoids those things. Inversion is replaced with solving a linear system, and the Jacobian is computed once before the loop. I also added argument-error checking and made TOL and itmax optional with default values. The solution is returned if convergence (as measured by TOL) is achieved; otherwise an error message is given.

restart
:

NewtonRaphsonSYS:= proc(
   X::list(And(name, Not(constant))),
   F::And(
         list(algebraic),
         satisfies(F->
            nops(F)=nops(X) and
            indets(F, And(name, Not(constant))) subset {X[]}
      )
   ),
   X0::And(list(complexcons), satisfies(X0-> nops(X0)=nops(X))),
   #optional keyword arguments and their default values:
   {  
      TOL::And(float,positive):= 10^(1+nops(X)-Digits),
      itmax::nonnegint:= 99
   }
)
local
   n, f:= <F>, Xn:= <X0>, H:= <seq(TOL, n= 1..nops(X))>,
   J:= VectorCalculus:-Jacobian(f,X)
;
   for n to itmax while LinearAlgebra:-Norm(H,2) >= TOL do
      H:= LinearAlgebra:-LinearSolve(
         rtable(evalf(eval(J, X=~ seq(Xn))), datatype= float),
         rtable(evalf(eval(f, X=~ seq(Xn))), datatype= float)
      );
      Xn:= Xn - H;
      userinfo(1, 'Newton', 'NoName', 'n'=n, 'Xn'=evalf[8]([seq(Xn)]));
      userinfo(1, 'Newton', 'NoName', `      `, 'H'=evalf[3]([seq(H)]))
   end do;
   if n > itmax then error "did not converge" end if;
   userinfo(1, 'Newton', 'NoName', `\n`); #blank line
   [seq(Xn)]
end proc
:   

#Usage example:
F:= [x*y-z^2-1, x*y*z+y^2-x^2-2, exp(x)+z-exp(y)-3];
X:= [x,y,z]:
X0:= [1,1,1]:

[x*y-z^2-1, x*y*z-x^2+y^2-2, exp(x)+z-exp(y)-3]

infolevel[Newton]:= 1: #Turn on userinfo statements.
Sol:= NewtonRaphsonSYS(X, F, X0, TOL= 1e-8, itmax= 10);

n = 1 Xn = [2.1893261 1.5984752 1.3939006]

       H = [-1.19 -.598 -.394]
n = 2 Xn = [1.8505896 1.4442514 1.2782240]
       H = [.339 .154 .116]
n = 3 Xn = [1.7801612 1.4244360 1.2392924]
       H = [0.704e-1 0.198e-1 0.389e-1]
n = 4 Xn = [1.7776747 1.4239609 1.2374738]
       H = [0.249e-2 0.475e-3 0.182e-2]
n = 5 Xn = [1.7776719 1.4239606 1.2374711]
       H = [0.279e-5 0.328e-6 0.270e-5]
n = 6 Xn = [1.7776719 1.4239606 1.2374711]
       H = [0.314e-11 0.422e-13 0.441e-11]

[HFloat(1.7776719180107403), HFloat(1.423960597888489), HFloat(1.2374711177317033)]

#Check residuals:
eval(F, X=~ Sol);

[HFloat(0.0), HFloat(-1.8388068845354155e-16), HFloat(-1.7763568394002505e-15)]

 


 

Download Newton.mw

The ​​​​​​solve command, used by Rouben, says that's there's no unconditional or "free" solutions, which is not suprising given that there are 12 equations for 6 unknowns. The eliminate command is an alternative that gives conditions under which solutions exist. I think that there's a good chance that that's what you're interested in. The following worksheet is meant to be put at the end of Rouben's. 
 

Sols:= eliminate(sys, vars):

These are the expressions that must be equal to 0 in order for the following solution to be valid:

Conds:= <Sols[2][]>;

Vector(6, {(1) = (p[1, 2]-p[1, 3])*alpha, (2) = -(p[2, 2]-p[2, 3])*alpha, (3) = p[1, 1]-p[1, 3], (4) = -p[2, 1]+p[2, 3], (5) = alpha*p[1, 2]-alpha*p[1, 3]-p[1, 1]+p[1, 3], (6) = -alpha*p[2, 2]+alpha*p[2, 3]+p[2, 1]-p[2, 3]})

(1)

And these are the solutions for the q's:

Sol:= < (<lhs~([Sols[1][]])> =~ ``) | <rhs~([Sols[1][]])>>;

"[[[q[1,1]=,-(beta p[2,3]-beta-p[1,3]+1)/((beta^2-1) (alpha+1))],[q[1,2]=,-(alpha beta p[2,2]-alpha beta p[2,3]-alpha p[1,2]+alpha p[1,3]+beta p[2,2]-beta-p[1,2]+1)/((beta^2-1) (alpha+1))],[q[1,3]=,-(beta p[2,3]-beta-p[1,3]+1)/((beta^2-1) (alpha+1))],[q[2,1]=,-(beta p[1,3]-beta-p[2,3]+1)/(alpha beta^2+beta^2-alpha-1)],[q[2,2]=,-(alpha beta p[1,2]-alpha beta p[1,3]-alpha p[2,2]+alpha p[2,3]+beta p[1,2]-beta-p[2,2]+1)/((beta^2-1) (alpha+1))],[q[2,3]=,-(beta p[1,3]-beta-p[2,3]+1)/(alpha beta^2+beta^2-alpha-1)]]]"

(2)

 


 

Download Conditional.mw

Numerous applications, particularly in number theory, require only the integer part of the base-2 logarithm. There is a special command for that, ilog2(6), which is much more efficient than floor(log[2](6)).

This is done with the identity option to solve.

solve(
   identity(
      cos(t)^6+a*cos(t)^4*sin(t)^2+b*cos(t)^2*sin(t)^4+c*sin(t)^6 = cos(6*t), 
      t
   ), 
   {a,b,c}
);
                   {a = -15, b = 15, c = -1}

 

These very small imaginary parts are due to round-off error in float computations. I call them spurious imaginary parts, because the true answers are real. They can be removed by

simplify(fnormal([%]), zero);

In some exceptional cases, you may need to specify a tolerance (a number of digits) as a second argument to fnormal. See help page ?fnormal.

The RealDomain package is not appropriate for this. It's better to generate the answers with the imaginary parts and then remove them than it is to try to generate them without the imaginary parts.

The command to turn an expression into a function is unapply.

Mat:= module()
option package;
uses CF= CurveFitting;
export 
     LinModel:= proc(xliste::list, yliste::list, $)
     local a, b, x;
          unapply(
               CF:-LeastSquares(xliste, yliste, x, 'curve'= a*x+b),
               x
          )
     end proc; # LinModel 
end module; #Mat

If you are returning a function, there's no need to pass x to LinModel. Thus, I've made it a local.

What you've called sys is already a set because you put it in braces; then you've put that and the initial condition into more braces. So, your first argument to dsolve isn't a set of equations; rather, it's a set that contains an equation and also contains a set of equations. So, you can make your dsolve command

dsolve({sys[], varphi(0) = 0});

The [] after sys extracts the elements of the set sys. I do not understand why you include the arguments Xa, Xb, Ya, Yb in the command.  

Substitute exact fractions such as beta = 3/10 rather than decimals such as 0.3. If you use a decimal before calling dsolve, it will automatically change it to an exact fraction, which is why your (4) is correct. 

L:= [a,b]:
`/`(L[]);

             a/b

If for some strange reason you want your convert command to work, do

`convert//`:= (L::list)-> `/`(L[]):

Now,
convert([a,b], `/`);
             a/b
works. I say "strange" reason because the first way shown seems so much easier than using convert. But in general, to make convert(..., foo) work, define a procedure named `convert/foo`; or if you're curious how a stock command such as convert(..., list) works, then look at showstat(`convert/list`).

 

It'll be a lot easier to display summation notation than "..." form, so how about

Sum(q[i], i= 1..n);

?

First 126 127 128 129 130 131 132 Last Page 128 of 395