Carl Love

Carl Love

28095 Reputation

25 Badges

13 years, 100 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

To get the slope, you extract the points from the plot, then you use linear least-squares regression on the logarithms of the points. Here's an extended example:

 

restart:

Digits:= 15:

f:= (t,y)-> I*y:

Exact:= unapply(rhs(dsolve({diff(y(t),t) = f(t,y(t)), y(0)= -I})), t);

(1)

Mid:= proc(h::{positive, realcons})
local
     y0:= -I, #Initial condition
     t,
     y1:= evalf(y0+h/2*(f(0,y0)+f(h, y0+h*f(0,y0)))), #one step of Heun
     y2
;
     for t from h by h to 1 do
          y2:= y0 + 2*h*f(t,y1);
          y0:= y1;
          y1:= y2
     end do;
     #Return the error
     evalf(abs(Exact(t-h)-y0))
end proc:
          

plots:-loglogplot(Mid, 1e-5..1e-1, labels= [h, `error`]);

 

Extract points data matrix from the plot.

XY:= op([1,1], %);

(2)

To get the slope of the line, perform linear least-squares regression on the logarithms of the points data.

Statistics:-LinearFit([1,x], ln~(XY), [x]);

(3)

So the slope is the coefficient of x.

 

Download Midpoint_error.mw

The unusual form of the solution is due to gamma's special role as a constant. Although I don't know the exact reason for the mysterious appearance of I in the solution, I do know that it can be avoided by declaring gamma as local:

restart;
local gamma;
gamma = 1/sqrt(1-beta^2);
solve(%, beta
);

It is very difficult for me to read the MapleMath "font" of this website. It would best if you didn't use it and just used plaintext instead. Too bad that we don't have MathJax (a LaTeX subset used on many websites). I think that what you've entered is the arclength formula:

int(sqrt(1+diff(f(x), x)^2), x= 0..17)

Is that correct?

I presume that your function f(x) is real valued. Is that correct? Then the imaginary part is spurious, perhaps the result of evaluating an elliptic function in the symbolic antiderivative.

You'll very likely get faster, more accurate, and entirely real results if you use numeric integration:

evalf(Int(sqrt(1+diff(f(x), x)^2), x= 0..17));

Note that Int is with an uppercase I for numeric integration.

If the total number of index values (i.e., mul(k_||j+1, j= 1..n)) is small enough that the Cartesian product can easily fit in memory (say, less than 10-100 million), then a non-iterated solution is probably much faster. Here's one:

CartProdSeq := proc(L::seq(list))
local Seq,i,j;
option `Copyright (C) 2007, Joseph Riel. All rights reserved.`;
     eval([subs(Seq= seq, foldl(Seq
                              , [cat(i,1..nargs)]
                              , seq(cat(i,j)=L[j],j=nargs..1,-1)
                             ))]);
end proc:

#K:= [k1, ..., kn]:
K:= [2,3,4]:
add(f(i[]), i in CartProdSeq(seq([$0..kj], kj in K)));

While the example seems silly to me, the output seems no more jumbled than any other polynomial in Maple. Polynomials are not generally sorted by powers in Maple.

Let expr be your arbitrarily complicated expression. Then

evalindets(expr, realcons, evalf[30]);

Change float to float[8].

Your code can be sped up significantly by using evalhf. It may be sufficiently faster to just replace evalf with evalhf. But it may be made even faster than that by filling the entire Array with four for loops inside a single call to evalhf.

It's a bit tricky because the number of entries in each row is not the same. We get around this by using an internal invocation of display so that the Vector constructor < > only sees one entry in the second row. Specifically, it can be done like this:

plots:-display(< H, plots:-display(< P | Q >) >);

The < H, plots:-display(< P | Q >) > specifies a column Vector whose first row is H and whose second row is the result of another display of a Vector. The < P | Q > specifies a row Vector whose entries are P and Q.

Here's the worksheet:


restart:

H:= plots:-animate(plot3d, [x-k*y+1, x=-10..10, y=-10..10], k=-10..0, frames=4):

P:= plots:-animate(plot, [sin(x+t), x=-Pi..Pi], t=-Pi..Pi, frames=8):

Q:= plots:-animate(plot, [cos(x+t), x=-Pi..Pi], t=-Pi..Pi, frames=8):

 

plots:-display(< H, plots:-display(< P | Q >) >);

 

Download array_animation.mw

 

Sorry, it seems that MaplePrimes will not allow me to upload an array of animations. But it is in the worksheet.

Here's a hack to do it with no mouse and no GUI. This involves using a custom operator &= instead of the usual :=.

restart:
`&=`:= proc(f::symbol, expr::uneval)
local x:= eval(expr);
     print(op(1, subs(_f= f, _x= x, proc(_f:= expr=_x) end proc)));
     assign(f,x)
end proc:

a:= 10:  b:= 3:
f&= (a+b^2);


There are two caveats to using this: The expression to the left of &= must be a simple symbol rather than an indexed name (f[1]) or a function (f(1)). And the expression to the right of &= must be in parentheses (unless it only has one part, such as a function call), because it's impossible to change an operator's precedence (see ?precedence ).

Let me know if that's satisfactory because I have a few other ideas if it's not.

See ?printlevel . The default value of the environment variable printlevel is 1. Each level of loop nesting (or if statement nesting) adds one level. So, if printlevel is set to 2 or more, then you will see what's happening inside your double loop.

Local variables and parameters evaluate differently from global variables. Quoting from ?eval (seventh paragraph of Description):

 The default evaluation rules in Maple are full evaluation for global variables, and one-level evaluation for local variables and parameters.

When your eq2 (in procedure f) is evaluated to one level, that leaves its unevaluated. So, the cure is to use eval(eq2), which forces a full recursive evaluation.

restart;
f:= proc()
local
     x,y,
     eq1:= 5+3*x=0,
     eq2:= 2+7*x-3*y-5*x*y=0
;
     x:= solve(eq1,x):
     y:= solve(eval(eq2),y):
lprint('x'=x,'y'=y);
end proc:

f();
x = -5/3, y = 29/16

To address your second question, about ||: The || operator simply was not intended to concatenate arbitrary expressions. It was intended as a convenient way to create indexed symbols. What you are trying to achieve can be easily done with the formatted print statement, printf:

printf("x is %a", -5/3);

The zip command is used to apply a two-parameter procedure to two matched arrays (or lists) of arguments.

restart:
B:= (x,y)-> fsolve(y = (1-exp(-x*b))/(1-exp(-50*b)), b):
X:= #Some array of x values
Y:= #Some array of corresponding y values
zip(B, X, Y);

Bisection is a method of last resort. Why do you want to use it? You can quickly get any number of roots of f(p) with RootFinding:-NextZero. After making Preben's unapply correction, do this:

N:= 100:  #Number of roots desired.
Root:= Array(1..N):
Root[1]:= -RootFinding:-NextZero(p-> f(-p), 0):
for n to N-1 do
     Root[n+1]:= -RootFinding:-NextZero(p-> f(-p), -Root[n])
end do:

The minus signs are because it seems that you want the roots in the negative direction.

All that you need is

series(y(x+h) - y(x-h), h=0, 3);

NumberToList:= proc(x::realcons)
uses S= StringTools;
     parse~(S:-Explode(S:-Select(S:-IsDigit, sprintf("%a", x))))
end proc:

First 317 318 319 320 321 322 323 Last Page 319 of 395