epostma

1429 Reputation

18 Badges

14 years, 117 days
Maplesoft

Social Networks and Content at Maplesoft.com

Maple Application Center
I am the manager of the Mathematical Software Group, working mostly on the Maple library. I have been working at Maplesoft since 2007, mostly on the Statistics and Units packages and on contract work with industry users. My background is in abstract algebra, in which I completed a PhD at Eindhoven University of Technology. During my studies I was always searching for interesting questions at the crossroads of math and computer science. When I came to Canada in 2007, I had nothing but a work permit, some contacts at Maplesoft (whom I had met at a computer algebra conference a year earlier), and a plan to travel around beautiful Canada for a few months. Since Maplesoft is a company that solves more interesting math and computer science related questions than just about anywhere else, I was quite eager to join them, and after about three months, I could start.

MaplePrimes Activity


These are replies submitted by epostma

@steweb : seeing your replies now (I had some things distracting me between starting to answer your question and finishing it), why don't you try this:

g:=proc(j)
local x;
  if j=0 then
    y->int(y*x,x=0..1)
  else
    y->int(y*g(j-1)(x),x=0..1)
  end if:
end proc:

Hope this helps,

Erik Postma
Maplesoft.

For those who want all the gory details, you can find them in the ?using_parameters help page, and also in these pages:

Hope this helps,

Erik Postma
Maplesoft.

For those who want all the gory details, you can find them in the ?using_parameters help page, and also in these pages:

Hope this helps,

Erik Postma
Maplesoft.

I would tweak your (nice) solution a little bit, Joe:

sys := { NULL
, diff(y(t),t) + y(t) = 0
, y(0) = 1
, s(0) = 0
, p(0) = 0.1
}:

integ := dsolve( sys
, 'numeric'
, 'events' = [[y(t)=p(t), [s(t)=t, p(t) = p(t)/10]]]
, 'discrete_variables' = [s(t), p(t)]
);

plots:-odeplot(integ, [t, p(t)], t=0..25, axis[2]=[mode=log]);
plots:-odeplot(integ, [t, s(t)], t=0..25);

This way p(t) holds the next "boundary" to be tested against, and s(t) holds the time at which the last boundary was achieved.

Hope this helps,

Erik Postma
Maplesoft.

I would tweak your (nice) solution a little bit, Joe:

sys := { NULL
, diff(y(t),t) + y(t) = 0
, y(0) = 1
, s(0) = 0
, p(0) = 0.1
}:

integ := dsolve( sys
, 'numeric'
, 'events' = [[y(t)=p(t), [s(t)=t, p(t) = p(t)/10]]]
, 'discrete_variables' = [s(t), p(t)]
);

plots:-odeplot(integ, [t, p(t)], t=0..25, axis[2]=[mode=log]);
plots:-odeplot(integ, [t, s(t)], t=0..25);

This way p(t) holds the next "boundary" to be tested against, and s(t) holds the time at which the last boundary was achieved.

Hope this helps,

Erik Postma
Maplesoft.

I like your solution, @Preben - one little tweak to make it a little more friendly for the user: you could include the evaluation call into f. In the first case, that can be done by defining f with

f := (k, y) -> eval(k * 's'(y));

and in the second case, using

f := (k, y) -> value(k * s(y));

In both cases, the kernel-initiated simplification kills the expression with the sum before the evaluation expands it.

If possible for your application, you could consider not including the value call at all at this stage and just work with the inert addition. Consider where the first point is that you absolutely need to evaluate the sum; for example, if by that time y has a known, floating point, value, then it will be much faster to evaluate the sum then than with symbolic y.

Hope this helps,

Erik Postma
Maplesoft.

I like your solution, @Preben - one little tweak to make it a little more friendly for the user: you could include the evaluation call into f. In the first case, that can be done by defining f with

f := (k, y) -> eval(k * 's'(y));

and in the second case, using

f := (k, y) -> value(k * s(y));

In both cases, the kernel-initiated simplification kills the expression with the sum before the evaluation expands it.

If possible for your application, you could consider not including the value call at all at this stage and just work with the inert addition. Consider where the first point is that you absolutely need to evaluate the sum; for example, if by that time y has a known, floating point, value, then it will be much faster to evaluate the sum then than with symbolic y.

Hope this helps,

Erik Postma
Maplesoft.

@prosolve : I think it's unlikely that there exists a closed form solution to that equation. However, for the application you suggest I think it would be perfectly fine to include a numerical solver into your integration algorithm. I converted the expression to rational numbers and replaced the functions with indexed variables to obtain

expr := 1686/396263*theta0[k+1]/DT^2+(-1/4*sin(theta0[k+1])*
lambda1[k+1]*DT^2-1686/396263*theta0[k]-1686/396263*DT*
omega0[k]-DT^2+1/4*cos(theta0[k+1])*lambda2[k+1]*DT^2)/DT^2

You could solve this in a few steps in a Newton solver, probably starting from the previous value of theta0[k+1] - which I assume is theta0[k]. Initially I thought I'd be fancy and try to get a better initial estimate than theta0[k], but I don't think it's worth it. (For what it's worth, what I did was:

expr2 := series(expr, theta0[k+1] = theta0[k], 2);
expr2 := eval(expr2, O = 0);
solve(expr2, theta0[k+1]);

which leads to

-(396263*sin(theta0[k])*lambda1[k+1]*DT^2+1585052*DT^2-396263*cos(theta0[k])*
lambda2[k+1]*DT^2+6744*DT*omega0[k]+6744*theta0[k]-396263*cos(theta0[k])*
lambda1[k+1]*DT^2*theta0[k]-396263*sin(theta0[k])*lambda2[k+1]*DT^2*theta0[k])
/(-6744+396263*cos(theta0[k])*lambda1[k+1]*DT^2+396263*sin(theta0[k])*lambda2[
k+1]*DT^2)

but just doing one Newton step is less work and will probably work better.)

Hope this helps,

Erik Postma
Maplesoft.

@prosolve : I think it's unlikely that there exists a closed form solution to that equation. However, for the application you suggest I think it would be perfectly fine to include a numerical solver into your integration algorithm. I converted the expression to rational numbers and replaced the functions with indexed variables to obtain

expr := 1686/396263*theta0[k+1]/DT^2+(-1/4*sin(theta0[k+1])*
lambda1[k+1]*DT^2-1686/396263*theta0[k]-1686/396263*DT*
omega0[k]-DT^2+1/4*cos(theta0[k+1])*lambda2[k+1]*DT^2)/DT^2

You could solve this in a few steps in a Newton solver, probably starting from the previous value of theta0[k+1] - which I assume is theta0[k]. Initially I thought I'd be fancy and try to get a better initial estimate than theta0[k], but I don't think it's worth it. (For what it's worth, what I did was:

expr2 := series(expr, theta0[k+1] = theta0[k], 2);
expr2 := eval(expr2, O = 0);
solve(expr2, theta0[k+1]);

which leads to

-(396263*sin(theta0[k])*lambda1[k+1]*DT^2+1585052*DT^2-396263*cos(theta0[k])*
lambda2[k+1]*DT^2+6744*DT*omega0[k]+6744*theta0[k]-396263*cos(theta0[k])*
lambda1[k+1]*DT^2*theta0[k]-396263*sin(theta0[k])*lambda2[k+1]*DT^2*theta0[k])
/(-6744+396263*cos(theta0[k])*lambda1[k+1]*DT^2+396263*sin(theta0[k])*lambda2[
k+1]*DT^2)

but just doing one Newton step is less work and will probably work better.)

Hope this helps,

Erik Postma
Maplesoft.

@Markiyan Hirnyk : Sorry, I must have missed the replies last month. If I have time I usually try to respond to such questions.

You asked how it compares to the exact solutions. You already gave most of the solution to that question, at least if you're mostly interested in eyeballing the result:

ss := dsolve(sys, f(x));
ss := solve(rhs(ss), x, AllSolutions);
ss := evalc(ss) assuming a > 1/4;
zz := (indets(ss, name) minus {a, constants})[1]:
p2 := plot(eval~(ss, zz =~ [seq(-4 .. -1)]), a = 1/4 .. 5, view=[DEFAULT, 0..15]);
plots:-display(p1, p2);

I assigned the plot in my original answer to p1; in the last line of the above code I show both that and the four corresponding exact solutions in one plot:

The ordering of the colours is the other way around here, so you can see that the correspondence is mostly quite good. At the far left there's a couple of points where the correspondence is less good; I haven't investigated it, but I suspect that that's mostly plotting inaccuracies and not inaccuracies of the computation method.

Hope this helps,

Erik Postma
Maplesoft.

@Markiyan Hirnyk : Sorry, I must have missed the replies last month. If I have time I usually try to respond to such questions.

You asked how it compares to the exact solutions. You already gave most of the solution to that question, at least if you're mostly interested in eyeballing the result:

ss := dsolve(sys, f(x));
ss := solve(rhs(ss), x, AllSolutions);
ss := evalc(ss) assuming a > 1/4;
zz := (indets(ss, name) minus {a, constants})[1]:
p2 := plot(eval~(ss, zz =~ [seq(-4 .. -1)]), a = 1/4 .. 5, view=[DEFAULT, 0..15]);
plots:-display(p1, p2);

I assigned the plot in my original answer to p1; in the last line of the above code I show both that and the four corresponding exact solutions in one plot:

The ordering of the colours is the other way around here, so you can see that the correspondence is mostly quite good. At the far left there's a couple of points where the correspondence is less good; I haven't investigated it, but I suspect that that's mostly plotting inaccuracies and not inaccuracies of the computation method.

Hope this helps,

Erik Postma
Maplesoft.

Hi MRI_user,

Glad that this works for you. Just a quick note - with this definition, MySecondFunction and MyMatrix are functionally equivalent. For any values u, v, w, we have that MySecondFunction(u, v, w) and MyMatrix(u, v, w) return identical matrices.

Erik.

Hi MRI_user,

Glad that this works for you. Just a quick note - with this definition, MySecondFunction and MyMatrix are functionally equivalent. For any values u, v, w, we have that MySecondFunction(u, v, w) and MyMatrix(u, v, w) return identical matrices.

Erik.

@Robert Israel : no, the idea would be to plot a formula that happens to render as the program plotting the formula, in the same spirit as Tupper's formula rendering as itself.

The item that makes Tupper's formula a bit of a cheat to me is that n, the number actually encoding the formula, is not part of what's displayed. With this formulation you'd fix that.

I guess one way of constructing the bitmap n in the program would be by using a font and a representation of the program, but I have a, probably unreasonable, hope that it might be possible to somehow construct that number through different means. But as you say, it's probably way too much work given the payoff.

Erik.

@Robert Israel : nice. In TTY maple, interestingly, some of the lines get "swallowed up" by other lines printed over top of them, but you can still read the lines that remain. In GUI it works beautifully of course.

I guess the more difficult challenge is to create a program that does the same without the use of (2D or 3D) TEXT plot structures...

Erik.

5 6 7 8 9 10 11 Last Page 7 of 21