acer

32385 Reputation

29 Badges

19 years, 336 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

It's not clear to me what you are trying to accomplish, since there is only one free variable in the example you gave.

Could you explain what you want to be the second independent variable? I am going to guess that you mean the upper value of the index of summation.

restart;

plot(add(BesselJ(n,r),n=0..5),r=-10..10);

plots:-contourplot(sum(BesselJ(n,r),n=0..i),
                   r=-10..10,i=1..5,grid=[201,5],
                   contourlabels=false,
                   filled, coloring=["Navy", "Orange"]);

 

 

Download contourplot_sum.mw

Or possibly this (though not so interesting...)

plots:-contourplot(sum(BesselJ(n,r),n=0..floor(i)),
                   r=-10..10,i=0..5,grid=[51,201],
                   contourlabels=false,
                   filled, coloring=["Navy", "Orange"]);

My oracle came up with these.

You'll have to judge whether the author's of these textbook examples might consider the solutions to be real, or what you think of squaring (or apply exp to)  both sides of an intermediate equation, etc.

restart;

#903
ode:=diff(y(x),x)-y(x) = x*y(x)^(1/2):
ic:=y(0)=4:
sol:=dsolve([ode,ic],y(x)):
H:=odetest(sol,ode):
simplify((rhs-lhs)(map(`^`,isolate(H,indets(H,radical)[1]),2)));
simplify(factor(subs(exp(x)=exp(x/2)^2,H))) assuming real;

0

0

restart;
#1966
#Book: Differential equations for engineers by Wei-Chau XIE, Cambridge Press 2010
#Section: Chapter 2
#Problem number: 92.
restart;
ode:=2*x*diff(y(x),x)-y(x) = ln(diff(y(x),x)):
sol:=dsolve(ode,y(x)):
H1:=odetest(sol[1],ode):
H2:=odetest(sol[2],ode):
simplify((rhs-lhs)(rationalize(map(exp,isolate(combine(H1),op(2,H1))))));
simplify((rhs-lhs)(rationalize(map(exp,isolate(combine(H2),op(2,H2))))));

0

0

restart;
#3035
#Book: Ordinary differential equations and their solutions. By George Moseley Murphy. 1960
#Section: Various 37
#Problem number: 1118.
restart;
ode:=(diff(y(x),x)^2+1)^(1/2)+a*diff(y(x),x) = x:
sol:=[dsolve(ode,y(x))]:
H1:=odetest(sol[1],ode):
H2:=odetest(sol[2],ode):
simplify((rhs-lhs)(map(`^`,isolate(H1,indets(H1,radical)[1]),2)));
simplify((rhs-lhs)(map(`^`,isolate(H2,indets(H2,radical)[1]),2)));

0

0

restart;
#3036
#Book: Ordinary differential equations and their solutions. By George Moseley Murphy. 1960
#Section: Various 37
#Problem number: 1119.
restart;
ode:=(diff(y(x),x)^2+1)^(1/2)+a*diff(y(x),x) = y(x):
sol:=[dsolve(ode,y(x))]:
H1:=odetest(sol[1],ode):
H2:=odetest(sol[2],ode):
simplify((rhs-lhs)(map(`^`,isolate(H1,indets(H1,radical)[1]),2)));
simplify((rhs-lhs)(map(`^`,isolate(H2,indets(H2,radical)[1]),2)));

0

0

restart;        
#3054
#Book: Ordinary differential equations and their solutions. By George Moseley Murphy. 1960
#Section: Various 37
#Problem number: 1146.
restart;
ode:=ln(diff(y(x),x))+4*x*diff(y(x),x)-2*y(x) = 0:
sol:=[dsolve(ode,y(x))]:
H1:=odetest(sol[1],ode):
H2:=odetest(sol[2],ode):
simplify((rhs-lhs)(rationalize(map(exp,isolate(combine(H1),op(2,combine(H1)))))));
simplify((rhs-lhs)(rationalize(map(exp,isolate(combine(H2),op(2,combine(H2)))))));

0

0

restart;
#3129
#Book: Differential Equations, By George Boole F.R.S. 1865
#Section: Chapter 7
#Problem number: 8.
restart;
ode:=x = a*diff(y(x),x)+(diff(y(x),x)^2+1)^(1/2):
sol:=[dsolve(ode,y(x))]:
H1:=odetest(sol[1],ode):
H2:=odetest(sol[2],ode):
simplify((rhs-lhs)(map(`^`,isolate(H1,indets(H1,radical)[1]),2)));
simplify((rhs-lhs)(map(`^`,isolate(H2,indets(H2,radical)[1]),2)));

0

0

 

Download examples_of_odes_that_do_not_odetest_ac_pt1.mw

[edit. 23/12/2020] Partially addressing a followup example. (End-points can be examined separately.)

# 3647
restart;
ode :=diff(y(x),x)=2*(x*sqrt(y(x))-1)*y(x):
ic  :=y(0)=1:
sol :=dsolve([ode,ic]):
res :=odetest(sol,ode):
solve(simplify(res)) assuming real;
                 RealRange(Open(-1), infinity)

Here is another way that happens to work for this example.

Using allvalues on the RootOfs is a better/stronger idea. It's a shame that a single solve call cannot find the exact real solution directly, even with its explicit option (and possibly with some mechanism to denote realness).

restart;

F := -1 - y - exp(x);

-1-y-exp(x)

G := x^2 + y*(exp(x) - 1);

x^2+y*(exp(x)-1)

Ty:=eliminate({F,G},{y});

[{y = -1-exp(x)}, {-(exp(x))^2+x^2+1}]

remove(has,map(s->solve(s union Ty[1],explicit),
              {solve({G} union Ty[2])}), RootOf);

{{x = 0, y = -2}}

Download solve_exp_poly_ex.mw

There are also some slightly more general but uglier variants which utilize frontend (as if temporarily freezing the exp call) and pass an inequality to solve to denote realness.

There is the iscont command. Its Help page is linked from the "See Also" part of the discont command's Help page.

There will always be problematic examples, however.

restart;

f := piecewise(x < 1, x, -x);

f := piecewise(x < 1, x, -x)

iscont(f, x=-infinity..1);

true

iscont(f, x=-infinity..1, closed);

false

iscont(f, x=1..infinity);

true

iscont(f, x=1..infinity, closed);

true

restart;

g := piecewise(x <= 1, x, -x);

g := piecewise(x <= 1, x, -x)

iscont(g, x=-infinity..1);

true

iscont(g, x=-infinity..1, closed);

true

iscont(g, x=1..infinity);

true

iscont(g, x=1..infinity, closed);

false

 

Download iscont.mw

There is also fdiscont, which can be finicky but may offer some recourse for examples for which symbolic analysis is problematic. It doesn't go well on your piecewise example, though.

The first thing is that you have a syntax mistake. You have,
    S*`union`(T intersect R);
which ought to be,
    S union (T intersect R);
It's possible that was an implicit multiplication mistake made in 2D Input mode.

restart;

lign1 := S union (T intersect R);

`union`(S, `intersect`(R, T))

lign2 := (S union T) intersect (S union R);

`intersect`(`union`(R, S), `union`(S, T))

L1, L2 := eval([lign1, lign2],[`union`=`&or`, `intersect`=`&and`])[];

`&or`(S, `&and`(R, T)), `&and`(`&or`(R, S), `&or`(S, T))

Logic:-Equivalent(L1, L2);

true

Download logic_set.mw

There are several ways to do this. Here is but one of them,

restart

V := RandomTools:-Generate(variable(length = 4))

yhcm

{`~`[parse](StringTools:-Explode(V))[]}

{c, h, m, y}

NULL

Download CAP13CombinatoireQuestion_ac.mw

I'm not sure why you want them in a set at the end (instead of a list or Vector). When you put them into a set you may lose duplicates.

You can remove all output (including plots) from a worksheet as follows:
  Maple 2015:  File->Edit->Remove Output->From Worksheet   (Alt-e  v  w)
  Maple 2020:  Evaluate->Remove Output from Worksheet  (Ctl-d)

That message about storing large data mostly refers to the output of large Arrays/MatricesVectors/rtables, when they are printed as output. I don't know how to reuse such stored results upon Close/Re-open. (It is stored in some text-encoded reference in the XML of the .mw file, potentially large.)

The preference for the large-computation-autosave-thingy is in a popup that apepars when you click the "Confirmation Dialogs" button on the Tools->Options->Interface popup tab. 

If you are using Explore then the embedded plot will lose the axis and viewing-range properties you put on it using the right-click menus, whenever the plot is subsequently updated.

Instead, you can use the view option in the actual plotting command to fix the visible range so that it remains static for all updates of the Sliders..

You may also use the sprintf command can construct strings. I think that often looks more legible than a complicated cat call.

You can escape double-quotes within the format. There are various ways to do it, depending on whether you want to pass allmembers as name or "allmembers" as string.

When you parse the result the inner double-quotes should survive. (Of course, parsing these results reuires Property to be a member of module NODESteelProfiles_I .)

g := "NODESteelProfiles_I:-Property";                                              
                        g := "NODESteelProfiles_I:-Property"

sprintf("%s(%a)", g, "allmembers");

             "NODESteelProfiles_I:-Property("allmembers")"

sprintf("%s(\"%a\")", g, allmembers);

             "NODESteelProfiles_I:-Property("allmembers")"

sprintf("%s(\"%s\")", g, "allmembers");

             "NODESteelProfiles_I:-Property("allmembers")"

The inner double-quotes are escaped within the string results. You don't see that when your merely print it, but you would if your used lprint.

The subsindets command does syntactic substitution. That is mostly similar to how subs works, and as the Help page explains by default subsindets will utilize subs to perform the replacement.

It is intentional that subsindets does not do a full evaluation by default. That allows multiple syntactic substitutions and structural manipulations to be performed in sequence.

You can compare subsindets with evalindets. You could also compare with calling the indexed command name subsindets[eval].

Entering an expression at the top-level does an evaluation (unless you use special syntax to control that). That's true even if done on something that subsindets or subs had returned unevaluated. It's also true if the expression is used for an assignment statement.

There should be nothing highly surprising in the following:

restart;

A := int(f(x), x);

int(f(x), x)

B := subsindets(A, specfunc(anything,f), u->op(u)^2);

int(x^2, x)

B;

(1/3)*x^3

eval(B, 1);

int(x^2, x)

evalindets(A, specfunc(anything,f), u->op(u)^2);

(1/3)*x^3

subsindets[eval](A, specfunc(anything,f), u->op(u)^2);

(1/3)*x^3

Download subsindets_evaluation.mw

As for int and your example involving exp(various ln calls, etc) it is much easier to snipe than to code. It's a difficult subject and there is no perfect algorithm for discerning when and whether to simplify all possible examples that might benefit as well as succeed in decent time. There is always room for improvement.

It seems to me that there are two kinds of issue worth considering here. One is the construction (binning, etc) of the data, and another is the plotting and visual rendering.

I'd like to look at the first part later, which I think could be done with some efficient procedure(s). There may be some related bits of functionality within the ImageTools:-Scale command. But I don't have time to study those R commands at the moment.

Here are a couple of alternatives for the plotting/rendering part, using techniques that are actually 2D.

You might still experiment with the nb size, naturally.

restart:

interface(version);

`Standard Worksheet Interface, Maple 2020.1, Linux, June 8 2020 Build ID 1474533`

with(Statistics):

UseHardwareFloats := false:

N := 10^4:
A := Sample(RandomVariable(BetaDistribution(2, 6)), N):
B := Sample(RandomVariable(BetaDistribution(6, 2)), N):

nb := 50:

mA := min(A):
hA := Range(A)/(nb-1):
mB := min(B):
hB := Range(B)/(nb-1):

cA := floor~((A-~mA)/~hA)+~1:
cB := floor~((B-~mB)/~hB)+~1:

NP := Matrix(nb$2, 0):
for n from 1 to N do
  NP[cA[n], cB[n]] := NP[cA[n], cB[n]]+1
end do:

NP:
 

# The OP's original, a 3D plot (be careful not to rotate it...)
# This is large, unweildy, and at large dimensions makes the GUI
# sluggish. In the worst case inadvertant mouse-selection may crash...

DPmat := plots:-matrixplot(
  NP,
  heights=histogram,
  colorscheme=["Blue", "White", "Red"],
  orientation=[270, 0, 0],
  style=patchnogrid,
  lightmodel=none,
  labels=["A", "B", ""],
  axis[1]=[tickmarks=[nb+1=max(B), 1=mB]],
  axis[2]=[tickmarks=[nb+1=max(A), 1=mA]]
):
DPmat;

# This produces a 2D structure akin to what the `densityplot` command
# returns.
# The values get interpolated (because the binned data is smaller than
# the dimensions at which we render here. That's not so nice because the
# effect is a blurred rendering.
# But this computation is a nice intermediate and very easy way to
# apply the colorscheme, and below we rip out the color Array from it.

DP:=plots:-surfdata(
  NP,
  mB..max(B), mA..max(A), dimension=2,
  colorscheme=["Blue", "White", "Red"],
  style=patchnogrid,
  labels=["A", "B"],
  axis[1]=[tickmarks=[mB, max(B)]],
  axis[2]=[tickmarks=[mA, max(A)]],
  size=[475,450]
):
DP;

RAW:=ImageTools:-Rotate(op([2],indets(DP,specfunc(anything,COLOR))[1]),'left'):

DPimg := plot(mB..max(B), mB..max(A),
     background=ImageTools:-Scale(RAW,1..400,method=nearest),
     size=[475,450], axes=boxed,
     labels=["A", "B"],
     axis[1]=[tickmarks=[mB, max(B)]],
     axis[2]=[tickmarks=[mA, max(A)]]):
DPimg;

length(DPmat), length(DP), length(DPimg);

1970649, 10408, 5425

ImageTools:-Embed(ImageTools:-Scale(RAW,1..400,method=nearest));

 

2D_density_plot_acc.mw

I'm not sure offhand why that happens. But you might convert the Matrix (or Vector...) to listlist before construction of the procedure.

restart;

with(CodeGeneration):

A := <<a|a^2>,<a^4|a^13>>;

Matrix(2, 2, {(1, 1) = a, (1, 2) = a^2, (2, 1) = a^4, (2, 2) = a^13})

P := unapply(convert(A,listlist), proc_options=[]);

proc () [[a, a^2], [a^4, a^13]] end proc

Matlab(P, resultname = "A");

function Preturn = P()
  Preturn = [a a ^ 2; a ^ 4 a ^ 13;];

Matlab(P, optimize, resultname = "A");

function Areturn = A()
  t1 = a ^ 2;
  t3 = t1 ^ 2;
  t5 = t3 ^ 2;
  Areturn = [a t1; t3 t5 * t3 * a;];

Download CG_Matrix_opt.mw

No, in your examples you do not have to declare the parameters of the inner procedure as locals of the parent procedure.

That is true whether you use arrow notation or proc...end proc. That is handled automatically. 

The situations in which it does matter include:

1) You use `unapply` instead, to create the procedure.

2) You are using `subs` to replace a parameter name within a procedure (thus used like a template). In this case I usually use a name beginning with double-underscore, since nobody should assign to that. This name is being replaced with something hopefully nicer, anyway.

ps. You showed two examples. The local `x` in the first example has no bearing on the `x` used as parameter to the arrow procedure. The second example is slightly better because it doesn't contain a wasted, unused local declaration. Otherwise the two act the same.

In you example you don't want assuming to be called until the replacement procedure (that forms the third argument passed to subsindets, which does the replacement) itself gets called.

You want assuming to be invoked within the body of the replacement procedure, not as part of the creation of that procedure.

Put it another way, you want the assuming called during the calls to the replacement procedure, rather than when that procedure is merely defined. When that procedure is defined the names present when it later gets called are not yet available.

One way in which this can be attained is by use of brackets in the code. I also show some other ways below.

restart;

expr:=(1+csgn(a)*a)/(3*csgn(b)*b);

(1/3)*(1+csgn(a)*a)/(csgn(b)*b)

The following doesn't work because of precedence.
(The failure is due to syntax.)

subsindets(expr,'specfunc( anything, csgn )', f->simplify(f) assuming positive);

(1/3)*(1+csgn(a)*a)/(csgn(b)*b)

Compare with the following.

subsindets(expr,'specfunc( anything, csgn )', f->(simplify(f) assuming positive));

(1/3)*(1+a)/b

The simplify call is not needed for this example.

subsindets(expr,'specfunc( anything, csgn )', f->(f assuming positive));

(1/3)*(1+a)/b

 

Here are some other variants.

 

subsindets(expr,'specfunc( anything, csgn )', ()->1);

(1/3)*(1+a)/b

subsindets(expr,'specfunc( anything, csgn )', f->`assuming`([f],[positive]));

(1/3)*(1+a)/b

subsindets(expr,'specfunc( anything, csgn )', f->`assuming`([f],[op(1,f)::positive]));

(1/3)*(1+a)/b

expr assuming positive;

(1/3)*(1+a)/b

indets(expr,'specfunc(csgn)');
map(`::`,op~(indets(expr,'specfunc(csgn)')),positive);
expr assuming map(`::`,op~(1,indets(expr,'specfunc(csgn)')),positive)[];

{csgn(a), csgn(b)}

{a::positive, b::positive}

(1/3)*(1+a)/b

Download subsindets_assuming.mw

I don't what all your other examples might be. But you may wish -- at some point -- to distinguish between the argument to csgn and the names in that argument (if it is not just a name).

There are several problems in the code. One of the main ones is that the procedures use the first argument f as if it were an operator/procedure, by applying it. But your example passes an expression instead of an operator/procedure.

Check for mistakes. (Vary the fixed step-size h, etc.)

The is the RK4 procedure, but you could edit the Euler version similarly.

The plotting below involves some recomputation (via RK4) that is not strictly necessary. I kept to what I perceived as the spirit of your original, and made no attempt to optimize plotting, etc. (I can think of at least four decent ways to accomplish better efficiency there: storing/returning all the iterates, memoization in a few ways, etc.)

restart

RK4 := proc (f, a, b, alpha, h, N) local t, w, i, k1, k2, k3, k4; t := a; w := alpha; for i to N do k1 := f(t, w); k2 := f(t+(1/2)*h, w+(1/2)*h*k1); k3 := f(t+(1/2)*h, w+(1/2)*h*k2); k4 := f(t+h, w+h*k3); w := evalf(w+(1/6)*h*(k1+2*k2+2*k3+k4)); t := t+h end do; w end proc

Digits := 15

15

F := unapply(exp(3-y), [t, y])

proc (t, y) options operator, arrow; exp(3-y) end proc

a, b, y0, h := 0, 3, .5, .1; N := floor((b-a)/h)

0, 3, .5, .1

30

RK4(F, a, b, y0, h, N)

4.12581955881883

RK4(F, a, b, y0, h, 0)

.5

exsol := simplify(dsolve({diff(y(t), t) = F(t, y(t)), y(a) = y0}, y(t)))

y(t) = 3+ln(t+exp(-5/2))

evalf(eval(exsol, t = b))

y(3) = 4.12560631556497

plots:-display(plot([seq([n*h, RK4(F, a, b, y0, h, n)], n = 0 .. N)], 0 .. 3, style = point, symbolsize = 15), plot(rhs(exsol), t = a .. b, color = blue))

 

Download rk4new_ac.mw

First 101 102 103 104 105 106 107 Last Page 103 of 336