Carl Love

Carl Love

16474 Reputation

23 Badges

7 years, 178 days
Mt Laurel, New Jersey, United States
My name was formerly Carl Devore. I was in the PhD math program at University of Delaware until 2005. I was very active in the Maple community at that time.

MaplePrimes Activity


These are answers submitted by Carl Love

There's no way that you're going to get an interesting result for this, because Maple doesn't even know how to differentiate KummerU with respect to its first parameter. But if you just want a result, any result, instead of an error message, you could do

series(' 'KummerU' '(p, 3/2, t), p);

That's two set of single quotes, not one set of double quotes.

Your underlying plot command is nonsense because it specifies a plot with respect to two variables, x1 and x2. There can only be one such variable.

Your trouble has nothing to do with the background option.

You can't specify both the value of the step size h and the desired accuracy at the same time because they depend on each other.

You'll also need to specify either the number of steps or the final x-value. Usually n is the number of steps, so are you sure that you're interpretting the instructions correctly?

Your dx and point_list don't make sense because they're loop invariants: Their value doesn't change in the loop.

I do this:

solve({sec(x)=sqrt(2), x>=3*Pi/2, x<=2*Pi}, allsolutions, explicit);

This'll also work with RealDomain:-solve. IMO, it's usually better to put inequalities inside the solve rather than in an assuming clause. (And the help ?RealDomain says that assumptions don't work with that package.)

This usage of explicit is undocumented, but I've been using it since Maple 16 (and it may have existed before that). I just discovered it by experimentation. Both the keywords allsolutions and explicit are required to get the specific results that satisfy the inequalities. Taking the help page ?solve,details at face value, it makes no sense to me that these keywords would work together for a transcendental equation, but the above results show that they obviously do. 

I would avoid using RealDomain unless you're sure that you need it. Wanting only real solutions is not a good reason to use it.

You can provide fsolve with ranges within which to search, if you know some:

fsolve({r1, r2}, {xa1= 0..1, xb1= 0..1});
      {xa1 = 0.5048762298, xb1 = 0.06850883394}

This solution is returned immediately.

Any system of recurrence relations of finite order with initial conditions (regardless of linearity, homogeneity, or autonomy) is equivalent to an ODE IVP system where the first derivatives are interpretted as the discrete derivative diff(x(n), n) = x(n+1) - x(n), and this is extended to higher-order derivatives in the natural way: (D@@k)(x)(n) = (D@@(k-1))(x)(n+1) - (D@@(k-1))(x)(n). The new equations are formally called difference equations rather than differential equations. If that ODE IVP is passed to dsolve(..., numeric, ...), and solved with Euler's method with stepsize = 1, then the solutions are identical to the solutions of the recurrence relations. These solutions can be plotted with plots:-odeplot, which makes this a very convenient method.

Here's a utility to convert a system of recurrence relations to its ODE IVP equivalent:

RecToDsys:= proc(
    Sys::set(algebraic=algebraic), #recurrence relations expressed with indexed names
    V::list(name), #dependent variables
    n::name, #independent variable
    ICs::set(indexed(integer)=complexcons) #initial conditions
)
local Z,K,k,j;
    (Z,K):= (min,max)(op~(lhs~(ICs)));
    evalindets(
        eval(Sys, n= n-min(eval(op~(indets(Sys, specindex(V))), n= 0))),
        specindex(V), 
        t-> 
            local f:= op(0,t), k, K:= eval(op(t), n= 0);
            add(binomial(K,k)*diff(f(n), [n$k]), k= 0..K)
    ) union 
    eval(
       `union`(
            seq(
                {(
                    (D@@k)~(V)(Z)=~ 
                    `~`[`+`](seq((-1)^(k-j)*binomial(k,j)*~index~(V,Z+j), j= 0..k))
                )[]},
                k= 0..K
            )
        ),
        ICs
    )    
end proc
:

In order to apply this to your system, I needed to make up 6 initial conditions. You could probably come up with better initial conditions than I could. 

Sys:= {
    x[n]=1/3*(2*x[n-1]*y[n-1]+4*x[n-1]*z[n-1])+1/12*(2*x[n-2]*y[n-2]+4*x[n-2]*z[n-2]),
    y[n]=1/3*(1/4*x[n-1]*z[n-1]+y[n-1])+1/12*(1/4*x[n-2]*z[n-2]+y[n-2]),
    z[n]=1/3*(x[n-1]*z[n-1]+2*y[n-1]*z[n-1])+1/12*(x[n-2]*z[n-2]+2*y[n-2]*z[n-2])
}:
Dsys:= RecToDsys(Sys, [x,y,z], n, {x[0]=-1, x[1]=0, y[0]=1, y[1]=0, z[0]=-1, z[1]=-2});
Dsol:= dsolve(Dsys, numeric, method= classical[foreuler], stepsize= 1):
Dsol(10);
plots:-odeplot(Dsol, [x,y,z](n), n= 0..10);

For the initial conditions that I tried, the solutions either quickly blew up (diverged to infinity) or quickly converged to 0. But you may know the proper initial conditions. My procedure does NOT require that the initial conditions start at index 0.

Simply replace a with a/((max-min)(a)/2). For example,

plot(<t | a/((max-min)(a)/2) >);

Use parentheses, not square brackets. 

Your equations are autonomous, which means that the independent variable n only appears as an index. And your equations are second order, meaning that each iteration can be computed strictly from the previous 2 iterations (and the independent variable if they weren't also autonomous). Just as with systems of differential equations, your system can be converted to an equivalent first-order system of 6 equations. At that point, you have an autonomous system of first-order recurrences, so this Post that I wrote 1-1/2 years ago should help you: "Autonomous first-order nonlinear recurrences with parameters and multiple dependent variables".

While you're studying that Post, I'll show how to convert your system to first order.

@Thomas Dean There should be a help page documenting the distributive properties of all operators over all other operators, and which happen automatically. There are two "partially invisble" operators:

  1. Function application: This can be thought of as a binary infix operator (without an apparent symbol) that takes a single left operand and an expression sequence (the contents of the parentheses) as its right operand. In prefix form, this operator is `?()`.
  2. Index application: Works the same as function application except that its right operand is the expression sequence in the [...]. Its prefix form is `?[]`.

And there are four "outfix" operators (i.e., the operators surround their operands):

  1. list constructor:  [...]: Its prefix form is `[]`. But if the [] have any left operand, it's index application, not list construction!
  2. set constructor: {...}: Its prefix form is `{}`.
  3. rtable row constructor: <...|...>: Its prefix form is `<|>`.
  4. rtable column constructor: <......>: Its prefix form is `<,>`.

And there is a "meta" postfix operator: the elementwise operator ~: Its prefix form is `~`[...], where the [...contains the prefix form of the operator that is the (left) operand of ~.

All of these are at least somewhat overloadable\overrideable, and that's quite commonly done for the rtable constructors.

I think that index application has the highest precedence of any operator. Function application is very high. If you've read ?evalapply, you now know that function application distributes over almost all other operators.

These help pages are very important for understanding operators:

  1. ?operators,precedence (although it doesn't cover the seven operators above)
  2. ?use
  3. ?evalapply
  4. ?object,operators

If you've tried a few solutions with dsolve, you've realized by now that it won't work without a specifying a numeric value of the parameter A. Ideally we'd like to learn how the character of the solution changes with different ranges for A. For this, use the parameters option to dsolve(..., numeric, ...):

restart:
ode:= diff(theta(t),t) = 1+A+(A-1)*cos(theta(t));
Q:= dsolve({ode, theta(0)= -Pi/2}, numeric, parameters= [A]):
PlotQ:= proc(A) Q(parameters= [A]); plots:-odeplot(Q, args[2..]) end proc
:
MultiPlot:= proc(As::list(realcons))
    local k, N:= nops(As)-1; 
    plots:-display(
        seq(PlotQ(As[k+1], t= -2..2, color= COLOR(HUE, .85*k/N)), k= 0..N), 
        _rest
    )
end proc
:
plots:-animate(
    MultiPlot, 
    [['seq'(k, k= K..K+2, .05)], title= 'sprintf'("A = %4.2f..%4.2f", K, K+2)], 
     K= -3..1, frames= 200, paraminfo= false
);

I made this a separate Answer because it's a completely different solution to the problem, which doesn't use textplot.

Starting exactly where your original worksheet left off

PrePlot:= %:
dsol:= dsolve(
    {sys11[], z(0)=z0, y(0)=y0}, 
    numeric, parameters= [z0,y0]
):
Dots:= proc(ics)
local k;
    dsol(parameters= eval([z(0),y(0)], ics));
    plot(
        [seq](eval([[z(t), y(t)]], dsol(k)), k= [-1,0,1]),
        color= [green, yellow, red], style= point,
        symbol= soliddiamond, symbolsize= 24, _rest
    )
end proc
:
plots:-display(
    [PrePlot, Dots(ics1, legend= [t=-1, t=0, t=1]), Dots(ics2)]
);

 

Every if statement must end with one of fiend, or end if. Every proc statement must end with either end proc or end. So, you can do this:

Fib:= proc(n)
option remember;
    if n=1 then return 1
    elif n=2 then return 2
    else Fib(n-1)+Fib(n-2);
    fi
end proc;

 

I'll answer the question for another "forked" recursion: computing binomial coefficients from the Pascal's Triangle formula,
C(n,k) = C(n-1, k-1) + C(n-1, k). The first example uses option remember; the second doesn't.

restart:
Bin:= proc(n::nonnegint, k::nonnegint)
option remember;
    if k > n then 0
    elif n=0 or n=k or k=0 then 1
    else thisproc(n-1,k-1) + thisproc(n-1,k)
    fi
end proc
:
CodeTools:-Usage(Bin(120,60));
memory used=0.75MiB, alloc change=0 bytes, cpu time=16.00ms, real time=7.00ms, gc time=0ns
              96614908840363322603893139521372656

restart:
#The only difference is that I commented out the "option remember":
Bin:= proc(n::nonnegint, k::nonnegint)
(* option remember; *)
    if k > n then 0
    elif n=0 or n=k or k=0 then 1
    else thisproc(n-1,k-1) + thisproc(n-1,k)
    fi
end proc
:
CodeTools:-Usage(Bin(25,12)); #Much smaller numbers than above
memory used=1.11GiB, alloc change=0 bytes, cpu time=10.56s, real time=10.71s, gc time=625.00ms
                            5200300

Now, you just need to do the same thing with your procedure Fib.

Note that without option remember, computation of the first 30 Fibonacci numbers will take about 2^30 (~1 billion) procedure calls. So, I recommend that you lower that from 30 to 20 or so.

The command plots:-textplot can be used to plot text at any coordinates, at any angle, and in any color. If that text were an appropriately placed ">", it would look like an arrowhead on the curve.

f:= x-> x^2:
(Xmin,Xmax):= (0, 3):
PrePlot:= plot(f(x), x= Xmin..Xmax, color= red):
(Ymin,Ymax):= (min,max)(op([1,1], PrePlot));
                      Ymin, Ymax := 0., 9.

Angle:= arctan@(D(f)/((Ymax-Ymin)/(Xmax-Xmin))):
plots:-display(
    PrePlot,
    plots:-textplot(
        [seq([k, f(k), ">", rotation= Angle(k)], k= ceil(Xmin)..floor(Xmax))],
        color= red
    )
);

For the angle to work correctly, Ymax-Ymin must be the length of the y-axis and Xmax-Xmin must be the length of the x-axis. The relationship of these quantities to the actual curve is irrelevant.

1 2 3 4 5 6 7 Last Page 1 of 285