Carl Love

Carl Love

28035 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

That is bizarre. I've never seen that __SELECTION, nor is it happening for me when I execute your code.

Here's a way around that. And, anyway, this is something you should do even if you don't get that weird output. Your function can be integrated once, for symbolic x, to obtain an algebraic expression in x, which can then be evaluated for any numeric x. Okay, that probably sounds complicated, but it's all done with one command: unapply.

g:= unapply(int(exp(-2*t^2), t= x-1..x) - int(exp(-2*t^2), t= x..x+1), x);
evalf(g(1));

If you know that you're always going to want to do evalf with g, then you might as well roll the evalf into g's definition, as in

g:= evalf@unapply(int(exp(-2*t^2), t= x-1..x) - int(exp(-2*t^2), t= x..x+1), x);
g(1);

By the way, it's obvious that g(0) = 0, so I hope that your purpose for doing this is something deeper than finding a root of g.

 

The freezing that you're experiencing from numtheory:-divisors(X) is due simply to Maple's attempt to format and display the lengthy output. It's not due to any computational inability. There are two things that you can do to prevent this:

  1. (strongly recommended) Go to menu Tools -> Options -> Precision -> "Limit expression length to" and put in a relatively low number like 100000. Apply globally.
  2. (not so great, because it's easy to forget) If you anticipate lengthy output, end the command with a colon.

Note that if the computation is done and you're stuck in the formatting stage or outputing stage, you're totally out of luck! This is much, much worse than being stuck in the computation phase (where you can kill only the problematic kernel through an OS program like Task Manager). To kill the display phase, you need to kill the GUI, which terminates your entire Maple session---all open worksheets---and any unsaved work in those worksheets will be lost.

It's incomprehensible to me how Maplesoft could let such a basic and demoralizing problem go unfixed year after year.

 

It's straightforward. You translate your symbols to Maple's VectorCalculus notation, subtract the right side from the left side, and simplify. If you get 0, it's proven. Like this (I've used G for your weird scalar function symbol):

restart:
with(VectorCalculus):
v:= x,y,z:
F:= VectorField(<f(v), g(v), h(v)>, cartesian[v]):
alias(G=G(v)):
Curl(G*F) - (G*Curl(F) + Gradient(G, [v]) &x F);
simplify(%);

 

I think that the error message is clear enough. The syntax A:-B searches module A's exports for something named B. On the other hand, the B without a module prefix is just accessing a local from a parent in the same manner as would happen in most programming languages, because exports are also local.

There are two exceptions, however, under which A:-B searches through all of A's locals, not just its exports: (1) If A is the keyword thismodule, or (2) if kernelopts(opaquemodules) is set to false.

The easiest programmatic solution would be to have the numbers in a plain text file, two numbers per line separated by a comma, and use ImportMatrix (see ?ImportMatrix). Many formats are supported, and you don't really need the commas, nor does the file really need to be plaintext.

I'm not sure what you mean by "then forgotten". Do you mean that the numbers must be forgotten (or erased, or deleted), as if they posed some sort of security risk? Or do you just mean that you don't care what happens to the numbers after you've processed them?

I was aware of the existence of Maple's goto because I've occasionally seen it (very rarely) in code that I've looked at with showstat. I'd advise you not to use it. It looks dangerous. Note that what you think of as the statement label, common_exit, is really just a name, and common_exit: is simply a statement. If you try to send something to a non-existent label, I wouldn't be surprised at all if the kernel crashed.

In my youth, I once crashed an entire Burroughs 7700 mainframe computer (supporting hundreds of simultaneous users) by using a goto to a nonexistent label in a compiled Algol program. I got chewed out big time by the management. Fortunately, this label misdirection was very difficult to accomplish because the compiler was very strict about it. For one thing, it required that labels be locally declared as labels.

I avoid long chains of if x=... elif x=... elif ... ... else ... fi because it bothers me that x needs to be re-evaluated at each level. Plus they look untidy to me. Rather, I use something that I call a dispatch table---a table whose entries are procedures and whose indices are the values of x, which need not be numbers. Example:

Dispatch:= table([
   10 = proc() print("x = 10"); ... end proc,
   12 = proc() print("x = 12"); ... end proc
]);
if assigned(Dispatch[x]) then
   val:= Dispatch[x]()
else
   print("Don't know what to do with x =", x);
   val:= undefined
fi;

Of course, I don't actually use print statements.



 

 

 

The easy way to make a neatly formatted table is to use a Matrix:

interface(rtablesize= 21): #Allow display of matrices up to size 21x21.
Matrix([seq([t, evalf(-sin(t)), -.5*sin(.5*t)], t = 1 .. 20)]);

That's all there is to it! Note that I changed your curly braces {} to square brackets []. Don't use curly braces in any situation where the order of the contents within the braces matters.

With a tiny bit more effort, you can add column headers:

Matrix([[`time`, `Wave A`, `Wave B`], seq([t, evalf(-sin(t)), -.5*sin(.5*t)], t = 1 .. 20)]);

(which explains why I made the rtablesize 21 instead of 20).

The above can also be done with a DataFrame (see ?DataFrame), but, in this case, I didn't think it worthwhile to go that route.

The vector field <x, y> is not the same thing as the scalar function x*y. What you need is a path integral, not a line integral:

VectorCalculus:-PathInt(x*y, [x,y]= Path(<4*t, 3*t>, t= 0..1));   
     20

Or you may prepend Student:- to that; the syntax is identical.

Okay, here's a draft of a worksheet that you can start using right now. Make sure to heed the comment (set in a box of exclamation points) about reducing the pixels if your display is less than 4k.

I refer to your Q(x) as forcingfunc. I set it to simply x in this line:

OneFrame:= Explore_BVP(forcingfunc= x, scenes= [[[x,w(x)]], [[x,u(x)]]]):

You change that x to any function that you want. You can also arrange and combine the plots any way that you want (currently limited to horizontally) using that scenes option that I made up. There can be any number of plot windows left to right. Each sublist is one plot window. Each sub-sublist is a curve. Each algebraic element is a plot axis. Each window can be plotted in 2 or 3 dimensions depending on the number of elements in the sub-sublist. So, for example, if you want w and w' on the left plot, u vs. w on the middle plot, and a phase portrait for u on the right plot, you'd do

scenes= [[[x,w(x)], [x, diff(w(x),x)]], [[w(x), u(x)]], [[u(x), diff(u(x),x)]]]

If the refresh rate of plots is too slow as you move the slider, or maybe if you're getting a lot of popup error messages, change numpoints= 10000 to a lower number.

Maple code author: Carl Love <carl.j.love@gmail.com>

4 June 2018

Using the Explore command to plot a parameterized BVP

restart:

sys_ode:=
   2*C__5*diff(w(x),x) - C__1*diff(u(x),x$2) - 2*C__4*diff(w(x),x$3) = Q,
   2*C__2*u(x) - C__1*diff(w(x),x$3) - 2*C__3*diff(u(x),x$2) = 0
:
ics:= w(0) = 0, D(u)(0) = 1, (D@@2)(w)(0) = 0, D(u)(100) = 1, (D@@2)(w)(100) = 2:
 

sys_bvp:= unapply({sys_ode, ics}, [Q, C__||(1..5)]):

Explore_BVP:= proc(
   {forcingfunc::algebraic:= 0,
    scenes::list(list(list)):= [[[x,w(x)]],[[x,u(x)]]]
   }
)
   proc(C__1, C__2, C__3, C__4, C__5)
   local f, sol;
      try
         sol:= dsolve(
            sys_bvp(forcingfunc, C__1, C__2, C__3, C__4, C__5), numeric,
            abserr= 1e2, maxmesh= 2^14
         );
         plots:-display(
            `<|>`(seq(plots:-odeplot(sol, f, numpoints= 10000), f= scenes)),
            #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            #!! Values on next line are measured in pixels! You must change    !!
            #!! these if you're not using a 4k display. The numbers are        !!
            #!! horizontal pixels and vertical pixels for the plots            !!
            #!! themselves (as a whole), not for the overall Explore window.   !!  
            #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            size= [3000,1500]  
         )
      catch:
         cat("Something went wrong for those parameter values.", lastexception)
      end try
   end proc
end proc:
        

OneFrame:= Explore_BVP(forcingfunc= x, scenes= [[[x,w(x)]], [[x,u(x)]]]):

Explore(
   OneFrame(C__1, C__2, C__3, C__4, C__5),

   C__1= 1.2e8..1.8e8,
   C__2= 2e9..10e9,
   C__3= 2e8..6e8,
   C__4= 0..5e7,
   C__5= 3e7..5e7,

   #The next line is relative, ("explore width"/"screen width")x100%,
   #so, there's no need to worry about pixels.
   widthmode= percentage, width= 100
);

 


 

Download Explore_BVP.mw

(This one is a FAQ (not that we have a list of such).)

The keyword itself (i.e., the part that's on the left side of = when it's passed) is an unassigned global name. Thus, it needs an unary :- prefix to disambiguate it from the parameter of the same name (as is always the case when a global and some other name are ostensibly the same: Prefix the global). So, you need to do:

procB :=proc(the_equation)
  print("inside procB, the equation is ", the_equation);
  procA(':-the_equation'=the_equation); 
end proc:

The unevaluation quotes are optional (although good safe coding practice), but the :- prefix is mandatory in ambiguous cases. Some people use the prefix all the time; I find that excessive and ugly. If you don't already know that you've used an ambiguous name in your own procedure, well, then you're in pretty bad shape.

Note that no number of unevaluation quotes can ever turn a parameter back into its parameter name, as this example shows:

proc(x) print('''x''') end proc(5);

This is because the substitution of passed arguments for parameters happens before the execution of any code in the procedure's body.

You need to put y=5 in parentheses, (y=5), to disambiguate.

None of Maple's relational operators =, <>, <, <=, subset, in, etc., are associative: They all require exactly two operands; no more, no less; one on the left, one on the right. So, the expression A = B = C is ambiguous: Does it mean (A=B) = C or A = (B=C)? It's up to you to specify.

You may argue: "Of course = is associative, as well as symmetric and transitive!" Sure, that's the way that it's generally used in mathematical writing. But in computer languages, including your beloved Mathematica, it tends to also be used antisymmetrically, for example, as a means of passing arguments. Even in mathematical writing, "Let x = a" means something different from "Let a = x".

You've accidentally discovered that unevaluation quotes can act like parentheses. Try not to rely on this (unless you're code golfing). The quotes are very ephemeral. Examine the resulting expressions with lprint, and you'll see the parentheses.

What is Q? Even if you know what Q is, it's clear from your output shown above that you haven't told Maple what Q is. Or were you hoping to get a solution for an arbitrary forcing function Q? [LMFAO] Maybe that's what Maple is trying to do....

I suspect that if a symbolic solution is possible at all, it'll be enormously complex with deeply nested integrals and RootOf expressions using your C__s as coefficients. And what's the point of having an enormously complex symbolic solution that's incomprehensible to the human mind? I'll admit that there is sometimes a point to that, but what's your reason? Chances are good that a parametric numeric solver, specific to this BVP system, would suit your purpose better than a symbolic solution. That numeric solver could be written in 10-15 minutes, including on-screen sliders to change the solution's plot as the parameters are varied.

 

The numerical solution of the original ODE system can be obtained easily, quickly, and straightforwardly:

restart:
F:= [S,E,J,R]:
IVP:= [
   diff(S(t),t) = mu*N - (beta*J(t)+mu)*S(t),
   diff(E(t),t) = beta*S(t)*J(t)-(mu+sigma)*E(t),
   diff(J(t),t) = sigma*E(t)-(mu+gamma1)*J(t),
   diff(R(t),t) = gamma1*J(t)-mu*R(t),
   S(0)=S0, E(0)=E0, J(0)=J0, R(0)=R0
]:
params:= [
   N= 5e7, beta=0.00001, mu=0.02, sigma=45.6, gamma1=73,
   S0=12500000, E0=50000, J0=30000, R0=37420000
]:
sol:= dsolve(IVP, numeric, parameters= lhs~(params), maxfun= -1):
sol(parameters= params):
plots:-display(
   <seq(
       plots:-odeplot(
          sol, [t,f(t)], t= 0..70, numpoints= 5000
       ),
       f= F
   )>, thickness= 3
);

But I don't know how to solve the system after you've applied "homotopy...".

In the code below, I used gamma1 for gamma and J for I because the symbols gamma and I are reserved in Maple for specific constants. If you really want to use gamma and I, it's easy to make a few adjustments to accomodate that.

restart:
Sys:= [
   S(k+1) = (mu*N - beta*Sum(S(l)*J(k-l), l= 0..k) - mu*S(k))/(k+1),
   E(k+1) = (beta*Sum(S(l)*J(k-l), l= 0..k) - (mu+sigma)*E(k))/(k+1),
   J(k+1) = (sigma*E(k) - (mu+gamma1)*J(k))/(k+1),
   R(k+1) = (gamma1*J(k) - mu*R(k))/(k+1)
]:
V:= [S,E,J,R]: #The unknown functions
ICs:= V(0) =~ [S0,E0,J0,R0]: #The symbolic initial values

K:= 100:        #highest value of k to compute and plot
params:= [
   N = 5e7,     #population
   mu = 0.02,     #?
   beta= 1e-5,     #?
   gamma1= 73, #?
   sigma= 45.6,  #?
   S0= 1.25e7,     #S(0)
   E0= 5e4,     #E(0)
   J0= 3e4,     #J(0)
   R0= 3.742e7     #R(0)
]:
assign(eval(ICs, params)); 
NumSys:= subs(Sum= add, unapply(eval(Sys, params), k)):
for _k from 0 to K-1 do assign(NumSys(_k)) od:
plot([seq([seq([k,f(k)], k= 0..K-1)], f= V)], legend= V, labels= [k, ``]);

Use this:

SP:= (n::posint)-> add(mul~(ifactors(n)[2])):
SP(12);

The [2] part remains the same regardless of n. It's an index into the list structure returned by ifactors. See ?ifactors, and note the (slight) distinction between it and the more commonly used ifactor. If you use ifactors  instead of ifactor, the format is more consistent, and you're less likely to encounter exceptional cases. The ifactor form is intended primarily for display, while the ifactors form is intended primarily for subsequent computation.

 

First 174 175 176 177 178 179 180 Last Page 176 of 395