Carl Love

Carl Love

28020 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

In a Reply to @nm above, I mentioned the need for the student to begin with "geometric intuition" for this general class of problem. Both the Answer by Kitonum and the one by acer show great use of this intuition. However, in cases where Green's Theorem can be applied, very little intuition is needed. The Theorem can be applied when the curves can be parameterized; then the area integral can be replaced by a line integral. In other words, only the boundary matters (which is rather amazing), no matter how reticulated it is. (It's as if I could measure the area of the pond in front of my house just by walking around it's boundary, never looking across the water).

The limits of integration are determined by the intersection points alone; there's no need to find roots of the cubic.

In the following transcription, I've deleted  much of the output consisting of very long algebraic numbers (some of very high degree). You can see it in the worksheet. The given solution is exact.

restart:
Digits:= 15:
eq1:= x^2/5 - y^2/2 = 1:  eq2:= y = x^3 - 7/2*x^2 + 2*x:

#Parameterizations (must use a different parameter for each):
C1:= table([x= sec(s)*sqrt(5), y= tan(s)*sqrt(2)]):
C2:= table([x= t, y= subs(x= t, rhs(eq2))]):

#Find intersection points:
solve({C1[x]=C2[x], C1[y]=C2[y]}, {s,t});

#Filter out any that are not real:
R:= select(type, {allvalues(%)}, set(name=realcons));
#Count the intersections:
nops(%);
                               2

#Great! Having only 2 intersections makes this especially easy!

#We don't need float values for the rest of this problem, but just out
#of curiousity...
Rf:= evalf(R)
    Rf := {{s = -0.540880221771561, t = 2.60840228710112}, 
      {s = 0.716707677716986, t = 2.96571533700056}}

#We can now plot directly rather than using implicitplot:
plot(
    [
        [C1[x], C1[y], s= eval(s,R[1])..eval(s,R[2])],
        [C2[x], C2[y], t= eval(t,R[1])..eval(t,R[2])]
    ],
    thickness= 3
);

#The plot shows a simple connected region. Good!

#The following integral will either give the area or its negative,
#I don't care which.
A:= 
    int(C1[x]*diff(C1[y],s), s= eval(s, R[1]) .. eval(s, R[2])) +
    int(C2[x]*diff(C2[y],t), t= eval(t, R[2]) .. eval(t, R[1]))
;
evalf(A);
                        0.75833670589040

#I think that A is in an extension field of the algebraic numbers generated by 
#including the ln of a few of them, but I'm not sure.

op~(0, indets(A, function));
                          {RootOf, ln}
#A is extremely big:
length(A);
                             52898

GreensThm.mw (Download worksheet)

P.S.: My plot above taken by itself doesn't exclude the possibility that one of the curves doesn't re-enter the region for some unshown values of or t, which would substantially complicate the problem. A broader plot is useful as acer suggests.

Please read carefully the section "Using Assumptions" in the help page ?solve,details.

Conclusion: Don't use assuming real with solve.

After what acer taught you yesterday, I'd've thought that you'd stop using all-variable assumptions such as assuming real for any command. 

This code uses the extremely fast package LinearAlgebra:-Modular. It is much faster than the compiled code, finishing the 2^25 case in under 2 minutes.

[Edit: The complete corrected code is in a Reply below. It does the 2^25 case in 0.2 to 0.3 seconds.]

BooleMobiusTransform:= proc(V::Vector(datatype= integer[8]))
uses LAM= LinearAlgebra:-Modular; 
local 
    nv:= numelems(V), n:= ilog2(nv), im:= Scale2(1, n-1), istep:= im, 
    jm:= 1, h:= Scale2(im, 1), istart, 
    R:= LAM:-Create(2, h, 0, integer[8]);
    LAM:-Copy(2, V, 1..-1, R, 1..nv);
    to n do 
        istart:= 1; 
        to jm do
            LAM:-AddMultiple(
                2, R, istart..im-1+istart, 
                R, istart+istep..im-1+istart+istep
            );   
            istart+= h
        od;
        im:= Scale2(im,-1); istep:= Scale2(istep,-1); 
        jm:= Scale2(jm,1); h:= Scale2(h,-1)
    od; 
    R
end proc
:
for n from 11 to 19 do
    V:= Vector(2^n, rand(0..1), datatype= integer[8]):
    lprint('n'=n);
    V1:= CodeTools:-Usage(BooleMobiusTransform(V))
od:
n = 11
memory used=0.58MiB, alloc change=0 bytes, 
cpu time=15.00ms, real time=9.00ms, gc time=0ns
n = 12
memory used=1.16MiB, alloc change=0 bytes, 
cpu time=16.00ms, real time=14.00ms, gc time=0ns
n = 13
memory used=2.35MiB, alloc change=0 bytes, 
cpu time=31.00ms, real time=21.00ms, gc time=0ns
n = 14
memory used=4.76MiB, alloc change=0 bytes, 
cpu time=47.00ms, real time=44.00ms, gc time=0ns
n = 15
memory used=9.63MiB, alloc change=0 bytes, 
cpu time=93.00ms, real time=85.00ms, gc time=0ns
n = 16
memory used=19.51MiB, alloc change=-1.00GiB,
cpu time=390.00ms, real time=358.00ms, gc time=93.75ms
n = 17
memory used=39.51MiB, alloc change=1.00MiB,
cpu time=437.00ms, real time=393.00ms, gc time=109.38ms
n = 18
memory used=80.01MiB, alloc change=2.00MiB, 
cpu time=797.00ms, real time=742.00ms, gc time=125.00ms
n = 19
memory used=162.01MiB, alloc change=3.00MiB, 
cpu time=1.62s, real time=1.50s, gc time=265.62ms

 

Often (but not always), Maple can symbolically integrate the Heaviside forms of piecewise expresssions that it can't do directly. Two branches of your proposed result needed small corrections (assuming I and Maple did everything correctly here).

restart:
interface(showassumed= 0):
assume(n::posint, j::nonnegint, j <= n, T > 0, alpha > 0);
h:= T/n:
psi:= t->
    local i:= op(procname);
    convert(
        if i=0 then piecewise(t<0, 0, t<h, 1-t/h)
        elif i=n then piecewise(t < T-h, 0, 1-(T-t)/h)
        else 
            piecewise(
                t < (i-1)*h, 0, 
                t < i*h,     t/h-(i-1), 
                t < (i+1)*h, (i+1)-t/h
            )
        fi,
        Heaviside
    )
:    
g:= proc()
local i,j,t;
    (i,j):= op(procname);
    simplify(
        int((j*h-t)^(alpha-1)*psi[i](t), t= 0..j*h)/GAMMA(alpha)
    )
end proc
:
g[0,j]() assuming additionally, j > 0;
             1         //       alpha            (alpha + 1)
    - ---------------- \\(j - 1)      (j - 1) + j           
      GAMMA(alpha + 2)                                      

                       alpha\  (-alpha)  alpha\
       + (-alpha - 1) j     / n         T     /


#Compare with your proposed result:
#Note that I added "-" in front of (j-1)!
simplify(
    % - 
    h^alpha/GAMMA(alpha+2)*
        (-(j-1)^(alpha+1) + j^alpha*(1-j+alpha))
);
                               0

g[j,j]() assuming additionally, j > 0;
                         (-alpha)  alpha
                        n         T     
                        ----------------
                        GAMMA(alpha + 2)

#Clearly, that's your proposed result.

g[i,j]() assuming i::posint, additionally, i > j, 
    additionally, i <= n;
                               0

g[i,j]() assuming i::posint, additionally, i < j;
       1         / alpha /           alpha            
---------------- \T      \(j - i - 1)      (i - j + 1)
GAMMA(alpha + 2)                                      

                              alpha            alpha            \ 
   + (-i + j + 1) (-i + j + 1)      + 2 (j - i)      (alpha + 1)/ 

   (-alpha)\
  n        /


#Note how I corrected your formula below:

simplify(
    % - 
    h^alpha/GAMMA(alpha+2)*
        ((j-i+1)^(alpha+1) + 2*(j-i)^alpha*(alpha+1) - (j-i-1)^(alpha+1))
);
                               0

 

It can be done like this (Warning: The key line of this code contains a large number of weird looking punctuation marks):

SpecialIndexPrint:= (e::uneval)->
    subsindets(
        convert(op(0, eval(e,1)), `local`)[eval([op](eval(e,1)))[]] = 
            eval(e),
        typeindex(symbol),
        n-> nprintf(
            cat(`#msub(mi(\"%a\"),mo(\"`, `%a ` $ nops(n)-1, `%a\"))`),
            op(0,n), op(n)
        )
    )
:
for i to 2 do  
    for j to 2 do 
        S[i, j] := 2*mu*varepsilon[i, j] - add((2*mu)/3*varepsilon[r, r]*delta[i, j], r = 1 .. 3); 
        print(SpecialIndexPrint(S[i, j])); 
    end do;
end do:

Yes, my procedure SpecialIndexPrint is more complicated than the Answers by Kitonum and Tom, but it does have a few advantages:

  • It can be used either inside or outside the loop.
  • It is not specific to the name S.
  • It works on any indexed symbols (with any number of indices).
  • It replaces the commas with spaces (although it'd be easy to add a correction for this to either of those other Answers).
  • It replaces only the commas that occur as index separators.

The prettyprinted results of the above cannot be adequately displayed here in plaintext form. You'll need to execute it in a Standard worksheet to see it. But it's exactly what you asked for.

To better understand how the above procedure works, use lprint instead of print, for example,

lprint(SpecialIndexPrint(S[1,2]));

`#msub(mi("S"),mo("1 2"))` = 2*mu*`#msub(mi("varepsilon"),mo("1 2"))`-2/3*mu*
`#msub(mi("varepsilon"),mo("1 1"))`*`#msub(mi("delta"),mo("1 2"))`-2/3*mu*
`#msub(mi("varepsilon"),mo("2 2"))`*`#msub(mi("delta"),mo("1 2"))`-2/3*mu*
`#msub(mi("varepsilon"),mo("3 3"))`*`#msub(mi("delta"),mo("1 2"))`

 

Maple doesn't "simplify" equations, and from what you've shown of Mathematica (in this and in previous Questions), it doesn't either. When you use simplify on an equation, it just simplifies the two sides separately. Thus, the following can be and should be used for all equations, regardless of whether they have square roots:

(numer@simplify@(rhs-lhs))(new_ode);
                       3                          
             / d      \        / d      \         
             |--- u(x)|  + 8 x |--- u(x)| - 8 u(x)
             \ dx     /        \ dx     /         

In American algebra education, the word simplify is not applied to equations; one simplifies expressions (and equations are not considered expressions), and one solves equations. I assume that this is the same in Kitonum's country, so I agree with what he said about you confusing simplify and solve.

However, note that my solution above does not use solve! So I think that it is a black-box solution such as you need. Instead, the equation is converted to an expression, which we thus assume is equal to 0. The expression is simplified, and since it equals 0, we're justified in further simplifying it to just its numerator.

You'll make your coding much simpler if you immediately process all input equations with numer@(lhs-rhs). Indeed, in this case, the simplify isn't even needed: (numer@(lhs-rhs))(new_ode) produces the same output. I just included the simplify for more general cases. I would guess that this process of converting equations into expressions is nearly universal in computer algebra. Certainly it's older than Maple.

Using the type &inside that I defined in the previous Answer, all you need is this:

evalindets(expr, specfunc(abs) &inside specfunc(ln), eval, abs= _@@0)

_@@0 is a fewest-characters representation of the multivariate identity function; you could also use abs= (x-> x).

Targetting only 1-argument abs is possible, but a little more complicated.
 

You can create your own structured types, and a generalization of what you describe is something that I want often enough that it's worth having a "prepositional" structured type for it. Specifically, I want to check for instances of type A inside type B.

TypeTools:-AddType(`&inside`, (e, `in`::type, out::type)-> e::out and hastype(e, `in`)):        
expr:=sin(x)+ln(abs(x))+ln(x+1/sqrt(abs(x+3)))+ln(x^3):

indets(expr, specfunc(abs) &inside specfunc(ln));

It's been several years since anything has been required as a first argument of specfunc. The first argument defaults to anything if none is supplied.

I do appreciate the ambiguities that acer pointed out; however, having read and Answered several hundred of your thoroughly written[*1] Questions (which I almost always Vote Up even if I don't Answer), I'm pretty sure that I understand what you mean:

  • You know the specifc meaning of function as a Maple type (and due to a recent Answer you know that the evaluated result of sqrt is not a function).
  • You want the recursive scan that indets (and its partner hastype) usually provide.

[*1] Although your English grammar, spelling, etc., could be improved, I still almost always understand what you mean. I do realize that English is probably not your native language. I'd be happy to give ad hoc grammar advice if that'd be appreciated.

Always use Matrix instead of matrix.

Please post your worksheet and state your Maple version.

I believe that you should use D, but only indirectly, called from your own procedure. Rather than loading D's remember table, my procedure uses a separate table that only holds your special derivatives. This handles both symbolic and numeric constants correctly. It's also a very short code. The fact that my table is declared sparse means that the derivatives of any symbols not specifically indexed (which are presumed to be constants) will be 0.

restart:
MyDerivatives:= table(sparse, [x= x*y, y= 1, z= z^2]):
MyD:= e-> evalindets(D(e), specfunc(D), d-> MyDerivatives[op(d)]):
  
MyD(x*y^2/2 + 3*y^3*z*b + 3 + a);
              1    3            2          3  2  
              - x y  + x y + 9 y  z b + 3 y  z  b
              2                                  

 

Here are the "under the hood" basics of how assuming works. Understanding this should clear up your confusion. It is obvious that the proper evaluation of F(G(e1, e2)) assuming a1, a2 requires that the assumptions be applied before F(G(e1, e2)) is evaluated. In other words, the normal evaluation rules (innermost to outermost) cannot be used. This command can be written in functional form as

`assuming`(F(G(e1, e2)), [a1, a2])

`assuming` is a procedure written in Maple; you can read its code with showstat. Although it's quite complicated code in general, you can see at the top that its first parameter is declared ::uneval

Returning to your original example, if for some reason you actually want to use the assumption only for the simplify and not for the odetest, and you want to do it in a single expression (i.e., without the intermediate assignment to res), use

simplify(odetest(mysol, ode) assuming []) assuming x>0

If you prefer, the [] could be replaced with nothing or anything or x::real or x::complex. The inner assumptions then supercede the outer when the inner expression is evaluated.

 

I think that the following will work in Maple 18. Some minor simplifications of the input syntax are available in more-recent Maple. My way of entering the system highlights that it is a matrix-vector linear system. There are other ways of entering this that others may consider simpler, although they do involve a bit more typing.

Is:= i__||(1..2);
sys:= {
    seq(diff([Is](t), t) =~ <1/2, -3; 2, -6>.<Is(t)> + <5*exp(-2*t), 0>),
    ([Is(0)] =~ [1,-1])[]
};
Sol:= dsolve(sys);

The exact solution that you get here looks more complicated than
it really is. That's because the coefficients are roots of
quadratic equations containing square roots. These numbers 
can be put in decimal form via:

evalf(%);

You also mentioned Maple Calculator. Many commands available in Maple are also available in Calculator. I suspect that dsolve is one of them, but I'm not sure. It's an extreme amount of work to enter alphabetic commands through the keyboard available in the Calculator. On the other hand, Calculator can be used to take a photo of the problem and translate it into Maple syntax which you can then upload to your Maple 18.

Name the two equations from (19) eq1 and eq2. Exactly as you've already entered them is fine; just begin the lines with eq1:= and eq2:= . There's no need to enter the next two equations; the solve command does that work:

Sol:= solve({eq1, eq2}, {U__mn, V__mn});

Then extract the coefficients like this:

ABs:= {
    seq(a[k] = coeff(eval(U__mn, Sol), W__mn, k), k= 0..2),
    seq(b[k] = coeff(eval(V__mn, Sol), W__mn, k), k= 0..2)
};

If you wish to assign the values of a[1], ..., b[3] (I'm not saying that it's a good idea to do this!!), do

assign(ABs);

 

 

When it's possible to solve a problem symbolically, it's usually useful to do that and then substitute the numeric values for the parameters. The symbolic solution provides invaluable insight into the physics of the problem. This problem can be solved symbolically by nearly any method that can be used numerically, including those posted by Kitonum and Tom. And here's another:

LinearAlgebra:-LinearSolve(
   Matrix(scan= band[1,1], [[280], rho*h*omega^2+~[100,444], [200]]),
   -C*<555, 6665>
);
                    [       /           2             \]      
                    [   5 C \111 h omega  rho - 217316/]      
                    [ - -------------------------------]      
                    [                 %1               ]      
                    [                                  ]      
                    [      /            2             \]      
                    [  5 C \1333 h omega  rho + 102220/]      
                    [- --------------------------------]      
                    [                 %1               ]      
                      2      4    2              2            
               %1 := h  omega  rho  + 544 h omega  rho - 11600

(A,B):= seq(eval(%, [h, C, rho, omega]=~ [3/100, 1, 7800, 222]));
                               -63994265     -549030931 
                     A, B := -------------, ------------
                             1330038150364  950027250260

evalf([%]);
                    [-0.00004811460858, -0.0005779107187]

 

Is there any indication in the documentation that algsubs might be appropriate for this?

Just use eval(ode, sol). Note that y(x) is necessarily a syntactic subunit[*1] of any expression in which it occurs, so there's no reason to think that algsubs might work better; and, as always with algsubs, there are many  reasons to think that it might work worse.

Why do so many people gravitate towards algsubs? and then after algsubs, they try applyrule? That those commands are houses built on sand can be guessed just by reading between the lines of their help pages. And if that doesn't convince you, read their code. The fundamental commands are subs and eval.

[*1] If you understand the representation of expressions as trees or directed acyclic graphs (DAGs), the syntactic subunits are the nodes of those graphs.

First 75 76 77 78 79 80 81 Last Page 77 of 395