## 1364 Reputation

12 years, 305 days
Maplesoft

## Social Networks and Content at Maplesoft.com

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.

## unapply...

For this case, the result can easily be obtained by evaluating the integral before the assignment to f: then indeed x is gone from f:

f := unapply(int(a*x, x=0..1), a);
int(f(x), x);

returns x^2/4. The difficulty comes when on the one hand you want to already perform the integration, but on the other hand you want to use a recursive definition that is evaluated on the fly.

Hope this helps,

Erik Postma
Maplesoft.

## Not sure what would cause g(0) to be une...

Hi steweb,

I'm not sure what would cause the phenomena you're describing; could you show us the full definition of g that you're using?

Re: your side note. I understand your point, but unfortunately we have to consider how the evaluation proceeds in order to get to a solution. You're correct that function application is not exactly the same as substitution, not in Maple and certainly not mathematically, but it's pretty close to what the Maple kernel does when a function call is evaluated. And if we want to use a computer to do math for us, there's not too much we can change about that. What we can change is which value is substituted in and what that value is substituted into.

If you define f := a -> int(a * x, x = 0 .. 1) then int is not yet called - the call to int is just recorded verbatim in f. Later, if you evaluate f(x), a already has the value x, so int receives the two arguments x^2 and x = 0 .. 1.

What you seem to have in mind is to first evaluate the integral and then evaluate the result at a = x. That process would be easy to do for simple expressions using unapply, but less so for more complicated situations where you want to do the definition recursively. So let's try to get absolutely clear what you want to do first.

Finally about your last question: making x a local in g ensures that there is nothing inside Maple that depends on that particular x, so in this case, once int gets called, its first argument is the product of the local x created for that particular invocation of g and something which cannot contain that x (unless you do your best to break it). x cannot be a local to int: you are passing it into int yourself, whereas local variables are created by Maple at invocation time and are available only from that point onwards.

Hope that helps,

Erik Postma
Maplesoft.

## Not sure what would cause g(0) to be une...

Hi steweb,

I'm not sure what would cause the phenomena you're describing; could you show us the full definition of g that you're using?

Re: your side note. I understand your point, but unfortunately we have to consider how the evaluation proceeds in order to get to a solution. You're correct that function application is not exactly the same as substitution, not in Maple and certainly not mathematically, but it's pretty close to what the Maple kernel does when a function call is evaluated. And if we want to use a computer to do math for us, there's not too much we can change about that. What we can change is which value is substituted in and what that value is substituted into.

If you define f := a -> int(a * x, x = 0 .. 1) then int is not yet called - the call to int is just recorded verbatim in f. Later, if you evaluate f(x), a already has the value x, so int receives the two arguments x^2 and x = 0 .. 1.

What you seem to have in mind is to first evaluate the integral and then evaluate the result at a = x. That process would be easy to do for simple expressions using unapply, but less so for more complicated situations where you want to do the definition recursively. So let's try to get absolutely clear what you want to do first.

Finally about your last question: making x a local in g ensures that there is nothing inside Maple that depends on that particular x, so in this case, once int gets called, its first argument is the product of the local x created for that particular invocation of g and something which cannot contain that x (unless you do your best to break it). x cannot be a local to int: you are passing it into int yourself, whereas local variables are created by Maple at invocation time and are available only from that point onwards.

Hope that helps,

Erik Postma
Maplesoft.

## Try this...

@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.

## Try this...

@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.

## More details...

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.

## More details...

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 determining the time at which this h...

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.

## For determining the time at which this h...

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.

## A slight tweak...

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.

## A slight tweak...

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.

## Include a numerical solver....

@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]);`

`-(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.

## Include a numerical solver....

@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]);`

`-(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.

## Mostly quite good...

@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.

## Mostly quite good...

@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.

 4 5 6 7 8 9 10 Last Page 6 of 20
﻿