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 replies submitted by Carl Love

@Shaaban 

KummerU(p, q, t) is (by definition) a solution of the ODE  t y'' + (q - ty' - p y = 0. The corresponding KummerM is the other solution. The independent variable of the ODE is t, so it's relatively easy to differentiate the solution, y, with respect to t

I said earlier that Maple didn't know how to differentiate it with respect to p. That's not 100% true. If you do FunctionAdvisor(KummerU(p,q,t)) and expand the section "differentiation rule", you'll see that there's a fantastically complicated formula of an infinite series of improper integrals and GAMMA functions given as the derivative with respect to the first parameter.

Do you have some reason to want the series with respect to p, or is it just curiosity?

It looks like you're trying to illustrate a linear programming problem. If so, what is the objective function, and do you want to minimize or maximize it?

Just to clarify that for other readers: The OP wants an example of using PDEtools:-PolynomialSolutions(..., HINT= ...).

Is the new larger file size a problem for you, or is it merely a curiosity?

@saketh I suspect that what you want isn't possible. Here's my reasoning by analogy: Suppose that you were computing a power-series solution to an ODE by, say, the Froebenius method. In all but some trivial cases, it would be impossible to compute the degree-n term without first computing all the lower degree terms.

@saketh Thank you. Your example completely illustrates your Question, which I now understand. Unfortunately, I have no answer, but maybe someone else does.

@9009134 I present here the two plots that we're talking about to make it easier for other readers who might be able to help with this. Both are constructed from 605 raw data points, and both have been "normalized" as we discussed above. 

This is "s1".

This is s2:

But could you define "dominant mode" for me? 

It appears that there is a horizontal shift of about 1/2 period between the two plots. Does this need to be taken into account?

Could you give a small example of a PDE which has two polynomial solutions, one of which is returned by PDEtools:-PolynomialSolutions and the other of which is the homogeneous polynomial that you want?

@notgoodatmaple This worksheet should help you understand many of the details from my (rather terse) Answer, although I haven't endeavored to explain every nuance of the syntax. Nonetheless, after you've experimented some with this worksheet, feel free to ask any specific questions about the syntax.
 

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..], gridlines= false)
end proc
:
MultiPlot:= proc(
    As::list(realcons), trange::range(realcons):= -2..2,
    {Legend::truefalse:= false}
)
    local k, N:= nops(As)-1;
    plots:-display(
        seq(
            PlotQ(
                As[k+1], t= trange,
                color= COLOR(HUE, .85*k/N),
                `if`(Legend, legend= [A=As[k+1]], NULL)               
            ),
            k= 0..N
        ),
        legendstyle= [location= right], _rest
    )
end proc
:

diff(theta(t), t) = 1+A+(A-1)*cos(theta(t))

The procedure PlotQ produces the plot for a single value of A, which is its first argument:

PlotQ(-2, t= -9..9, thickness= 3);

You can interogate what the current parameter value is:

Q(parameters);

[A = -2.]

As as there is a numeric value for the parameter, you can request specific numeric values of theta(t):

Q(-4), Q(4), Q(-9), Q(9), Q(99), Q(-99);

[t = -4., theta(t) = HFloat(1.9104993102616643)], [t = 4., theta(t) = HFloat(-1.9106294563198518)], [t = -9., theta(t) = HFloat(1.9106329863105198)], [t = 9., theta(t) = HFloat(-1.9106337351769633)], [t = 99., theta(t) = HFloat(-1.9106329949645597)], [t = -99., theta(t) = HFloat(1.910633435519435)]

PlotQ(0, t= -5..5, thickness= 3, color= red);

Oviously, A=1 is a straight line:

PlotQ(1, t= -5..5, thickness= 3, color= orange);

PlotQ(2, t= -5..5, thickness= 3, color= yellow);

PlotQ(4, t= -5..5, thickness= 3, color= green);

The procedure MultiPlot plots several curves on the same axes, using several specified values of A. The curves are color coded using the visible spectrum (red, orange, yellow, green, blue, violet and all continuous gradations between), so the first specified A is red, the last is violet, and the others are somewhere in between.

MultiPlot([-2, -1, 0, 1, 2, 4], -5..5, thickness= 3, Legend);

The following does the same thing for -2 <= A <= 4, spaced by 0.05

MultiPlot([seq(-2..4, 0.05)], thickness= 0);

To produce the animation shown in my Answer, I used 200 frames, each frame being a call to MultiPlot with 40 A values evenly spaced over an interval of width 2. Here's a smaller animation:

plots:-animate(
    MultiPlot,
    [['seq'(A_red..A_red+1, .2)], thickness= 3],
    A_red= -3..3, frames= 60
);

 


 

Download ODE_parameter_plot.mw

 

Considering how often Questions about fsolve returning unevaluated appear here on MaplePrimes, it would be appropriate if it returned a message akin to the "solutions may have been lost" that solve returns. Then the help page could advise on things to try when that message appears. This change can be made while still maintaining backwards compatibility.

@mehdi jafari If you have initial conditions that make this system interesting, I'd be curious to see them.

@Carl Love The bug is fixed. Either download the code again, or simply change the exponent of -1 from j to k-j.

@mehdi jafari Thank you. There's a small bug in that code regarding the alternation of signs in computing the initial conditions. It'll probably take me a few minutes to fix.

By "normalize", do you mean that you want to rescale the amplitude to [-1, 1]?

Hint: What's equivalent to the remainder of p(x)/(x-a) for a polynomial p?

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