Joe Riel

9530 Reputation

23 Badges

20 years, 28 days

MaplePrimes Activity


These are answers submitted by Joe Riel

See the help page ?worksheet,documenting,styles.

Not a proof, but may lead you to one

convert(tan(phi)*tan(phi+Pi/2), 'sincos');
                                                                      -1

There is no good reason to recompute D(f) and (D@@2)(f) inside the loop.  Better to assign them to local variables.  You may then notice that there is a problem, no argument is passed to (D@@2)(f).  If you also assign some intermediate computations, you can then quickly see where things are going wrong by tracing the procedure call (i.e. call debug(NEWTON_RALPH) before evaluating.
 

What is is 'n' doing there?  It is in the snippet ( f(x)^2) n

 

z := 0.7*x + 5e-10*y; 
fnormal(z);                                
                                                                    0.7 x

Here's one approach,

proc()
local dsys, data, w;
    dsys := {diff(v(t),t) = -x(t), diff(x(t),t)=v(t), x(0)=1, v(0)=0};
    data := dsolve(dsys, 'numeric', 'output'=Array([seq(0..2, 0.1)]));
    w := 15;
    printf("%*s\n", w,  map(convert, data[1,1], string));
    printf("%*f\n", w, data[2,1]);
end proc():
              t            v(t)            x(t)
       0.000000        0.000000        1.000000
       0.100000       -0.099833        0.995004
       0.200000       -0.198669        0.980067
       0.300000       -0.295520        0.955336
       0.400000       -0.389418        0.921061
       0.500000       -0.479426        0.877583
       0.600000       -0.564643        0.825336
       0.700000       -0.644218        0.764842
       0.800000       -0.717356        0.696707
       0.900000       -0.783327        0.621610
       1.000000       -0.841471        0.540302
       1.100000       -0.891207        0.453596
       1.200000       -0.932039        0.362358
       1.300000       -0.963558        0.267499
       1.400000       -0.985450        0.169967
       1.500000       -0.997495        0.070737
       1.600000       -0.999574       -0.029200
       1.700000       -0.991665       -0.128844
       1.800000       -0.973848       -0.227202
       1.900000       -0.946300       -0.323290
       2.000000       -0.909298       -0.416147

As a bous, with Maple13 you can use convert~(data[1,1], string) to convert the elements in the Vector stored in data[1,1] to a string.  A slightly more interesting method is to assign an inline operator for convert and then use it with the ~ operator:

`&c` := convert:
data[1,1] &c~ string;

You might read www.mapleprimes.com/blog/joe-riel/maple-overload, which discusses some of the issues with overload and how to best handle them.  If you do use overloading, I suggest you use the following general format (here given as an example module):

mymodule := module()
export g;
local f, InvalidInput;
    InvalidInput := proc() cat("\n  ", StringTools:-FormatMessage(lastexception[2..-1])) end proc;
    g := overload ( [NULL
                     , proc(R::list , S::list, $)
                       option overload;
                           R, S;  # verify all parameters
                           try
                               print("this is the version with TWO lists");
                               (* actual code *)
                               f(1/2);  # an internal call with invalid input
                           catch "invalid input:" :
                               error InvalidInput();
                           end try
                       end proc

                     , proc(R::list, $)
                       option overload;
                           R; # verify all parameters
                           try
                               print("this is the version with ONE list");
                               (* actual code *)
                           catch "invalid input:" :
                               error InvalidInput();
                           end try
                       end proc
                    ]);
    f := proc(x::integer) x/3 end proc:
end module:

mymodule:-g([],[]);
                     "this is the version with TWO lists"

Error, (in g) 
  invalid input: f expects its 1st argument, x, to be of type
integer, but received 1/2

The reason for the try/catch statements, with the call to InvalidInput that slightly modifies the error message, is to trap and reraise error messages of the form "invalid input:" but in a form that does not trigger the overload switch mechanism.  To see this, comment out the call to InvalidInput in the first procedure; you will then get

mymodule:-g([],[]);
                     "this is the version with TWO lists"

Error, invalid input: no implementation of g matches the arguments in call,
g([],[])

 

 

 

If you are using Maple input, and the code is all in one input block, then you can use the multi-line Maple comment, (* ... *)

> (*
>    These lines are commented out.
> *)

 

This won't be the easiest approach, but it may do the job.  I've arbitrarily assigned expressions for f, g1, and g2 in terms of x, y1, and y2.

alias(x=x(y1,y2)):

f := x^2 + y1/y2:
g1 := x*y1:
g2 := x+y2:

eq1 := f = 0:

d1 := solve(diff(eq1,y1), {diff(x,y1)});
                                   d          1
                           d1 := {--- x = - ------}
                                  dy1       2 x y2

d2 := solve(diff(eq1,y2), {diff(x,y2)});
                                    d        y1
                            d2 := {--- x = -------}
                                   dy2           2
                                           2 x y2

d3 := solve(diff(g1,y1)=0, {diff(x,y1)});
                                    d         x
                            d3 := {--- x = - ----}
                                   dy1        y1


d4 := diff(subs(d1,d3),y2);
                                 d                   d
                                --- x               --- x
                                dy2        1        dy2
                     d4 := {1/2 ----- + ------- = - -----}
                                 2            2      y1
                                x  y2   2 x y2

d5 := frontend(solve, [d4, diff(x,y2)], [{Not(specfunc(anything,diff))},{}]);
                              d              x y1
                      d5 := {--- x = - -----------------}
                             dy2              2
                                       y2 (2 x  y2 + y1)

d6 := {diff(g2,y2)=0};
                                   / d   \
                            d6 := {|--- x| + 1 = 0}
                                   \dy2  /


d7 := subs(d5, d6);
                                     x y1
                      d7 := {- ----------------- + 1 = 0}
                                      2
                               y2 (2 x  y2 + y1)

frontend(solve, [d7,y2], [{Not(identical(x))},{}]);
                      2      3    1/2                   2      3    1/2
             -y1 + (y1  + 8 x  y1)              y1 + (y1  + 8 x  y1)
       {y2 = ------------------------}, {y2 = - -----------------------}
                          2                                 2
                       4 x                               4 x

Alternatively, you can use declare the type of the parameter.  While that may not always be possible, it is generally more convenient to let Maple do the parameter processing.  For example

f := proc(x::integer) x^2; end proc:
f(1/2);
Error, invalid input: f expects its 1st argument, x, to be of type integer, but received 1/2

Here's a simple helper procedure that implements Robert's suggestion

p1 := plot(sin):
p2 := plot(cos):

recolor := proc(p, color)
    subsindets(p
               , 'specfunc(anything, {COLOR, COLOUR})'
               , clr -> COLOUR(RGB,op(ColorTools:-NameToRGB(convert(color,string))/255.0))
              )
end proc:

plots:-display(recolor~([p1,p2],[red,blue]));

Did you try the help page?  As an example,
 

taylor(x^3+1, x=2);
                                              2          3
                    9 + 12 (x - 2) + 6 (x - 2)  + (x - 2

Followup:  taylor returns a Maple type series.   Use convert/polynom to convert back to a polynomial.

convert(%,polynom);
                                            2          3
                      -15 + 12 x + 6 (x - 2)  + (x - 2)

Note that Maple automatically expands the linear term.

expand(%);
                                     3
                                    x  + 1

 

Jakubi's approach is simpler, but this may be of interest.  The basic idea is to use simplify/polar; alas, it does not currently handle assumptions so we will surgically repair it (temporarily).  Here's the limitation

simplify(polar(r,theta)^k) assuming k::real;
                                           polar(r,theta)^k

Here's the minor surgery

`simplify/polar` := subs('{numeric, realcons}' = 'satisfies'(rcurry(is,real))
                         , eval(`simplify/polar`)
                        ):

Now

simplify(polar(r,theta)^k) assuming k::real;
                                                  polar(r^k,theta*k)

With that fixed, we can now convert to polar form and use the new simplification.

zr := a+b*I;
zp := polar(zr);
zc := conjugate(zp);
y1 := zp^k + zc^k;
simplify(y1,polar) assuming real;
evalc(%);
                 2*((a^2+b^2)^(1/2))^k*cos(arctan(b,a)*k)

With Maple 13 you can append the tilde to a symbol and it will map over elements.  Thus

fg~(L);

would also work.  See ?operators/elementwise.

What do you expect to get?  Does p actually represent a general expression in a(t) and b(t)?

Consider

p := (a,b) -> a^2 + sin(b):
y := p(a(t),b(t)):

The frontend procedure can be used here to compute the derivative with respect to a(t):
 

frontend(diff, [y,a(t)]);
                                                    2 a(t)

Actually, that only worked because the sin didn't contain a(t). Consider
 

y := a(t) + b(t)*sin(a(t)+b(t)):

To use frontend effectively here you need to use its optional third argument that specifies what not to freeze.

frontend(diff, [y, a(t)], [{`*`,`+`, And(function,Not(identical(a(t),b(t))))},{}]);
                                          1 + b(t) cos(a(t) + b(t))

Suppose that your expression contains an inert function of a(t) and b(t).

frontend(diff, [f(a(t),b(t)), a(t)], [{`*`,`+`, And(function,Not(identical(a(t),b(t))))},{}]);
                                          diff(f(a(t), b(t)), a(t))

%;
Error, invalid input: diff received a(t), which is not valid for its 2nd argument

 

That raises an error because the diff procedure cannot deal with a function as the second argument. The proper way to express the derivative is with the D operator.  That can be done by using  convert/D inside the call to frontend,

frontend(`convert/D`@diff, [f(a(t),b(t)), a(t)], [{`*`,`+`, And(function,Not(identical(a(t),b(t))))},{}]);
                 D[1](f)(a(t), b(t))

 

First 79 80 81 82 83 84 85 Last Page 81 of 114