I know that your issue has already been solved by the other replies, but just for completeness, here is a possible solution involving the taylor command. As you already found out, the issue is that what you've written doesn't distinguish between the general variable, x, and the specific value for it, 3, once Maple enters the procedure f that you've written. One way to solve that is by using a variable instead of the value when you call f:

(**) expr := f(a);
                                            2    3
                          expr :=  1 + a + a  + a

(**) eval(%, a = 3);
                                      40

This gives you access to the formal expression, expr, first, which you can then evaluate at a point. If all you're interested in is the evaluation, you could use a local variable in f; then, you'll need to write it explicitly as a procedure instead of using arrow notation:

newF := proc(a)
local x;
  return eval(convert(taylor(1/(1-x), x=0, 4), polynom), x = a);
end proc;

This way, you know that x will just be a variable (with no value assigned to it) when taylor gets to do its thing. Now newF(3) will return 40 directly. (And actually, the convert/polynom step is unnecessary now, because evaluating a series at a point drops the big-Oh term.) Note that newF(a) still returns the formal expression.

Hope this helps,

Erik Postma
Maplesoft.

I know that your issue has already been solved by the other replies, but just for completeness, here is a possible solution involving the taylor command. As you already found out, the issue is that what you've written doesn't distinguish between the general variable, x, and the specific value for it, 3, once Maple enters the procedure f that you've written. One way to solve that is by using a variable instead of the value when you call f:

(**) expr := f(a);
                                            2    3
                          expr :=  1 + a + a  + a

(**) eval(%, a = 3);
                                      40

This gives you access to the formal expression, expr, first, which you can then evaluate at a point. If all you're interested in is the evaluation, you could use a local variable in f; then, you'll need to write it explicitly as a procedure instead of using arrow notation:

newF := proc(a)
local x;
  return eval(convert(taylor(1/(1-x), x=0, 4), polynom), x = a);
end proc;

This way, you know that x will just be a variable (with no value assigned to it) when taylor gets to do its thing. Now newF(3) will return 40 directly. (And actually, the convert/polynom step is unnecessary now, because evaluating a series at a point drops the big-Oh term.) Note that newF(a) still returns the formal expression.

Hope this helps,

Erik Postma
Maplesoft.

Hi,

It looks like you're better off making dTvalues a list instead of a vector. Replace your definition by

dTvalues := [seq(seq([d, T], d=-30..30, 10), T=1..20)];

and then you should use nops instead of rtable_dims (since dTvalues is not an rtable anymore), so that would look like

i=1..nops(dTvalues)

where you had

i=rtable_dims(dTvalues)

before.

The other worksheet looks like quite a different problem and I'm not sure I have the time now to look at it. Since it's a different problem, you could try posting a new post for that and seeing if someone else can help you with that.

Hope this helps,

Erik Postma
Maplesoft.

Hi,

It looks like you're better off making dTvalues a list instead of a vector. Replace your definition by

dTvalues := [seq(seq([d, T], d=-30..30, 10), T=1..20)];

and then you should use nops instead of rtable_dims (since dTvalues is not an rtable anymore), so that would look like

i=1..nops(dTvalues)

where you had

i=rtable_dims(dTvalues)

before.

The other worksheet looks like quite a different problem and I'm not sure I have the time now to look at it. Since it's a different problem, you could try posting a new post for that and seeing if someone else can help you with that.

Hope this helps,

Erik Postma
Maplesoft.

Or, as a slight improvement, in this case you can also use ?LinearFit , since the expression is linear in the parameters (p[1] and p[2]). That will generally lead to a more accurate result and use less time. In fact, if you use ?Statistics[Fit] , then Maple will decide whether LinearFit is applicable or not; you could call it like this, for example:

independentDataMatrix := <<IC50> | <REDOX>>;
H := Statistics[Fit](p1 * x1 + p2 * x2 + p3, independentDataMatrix, LIPO, [x1, x2]);

You can add the output=solutionmodule option, like HerClau did above, if you like. Otherwise it will just give you the fitted function.

You need to make a decision which quantities you consider dependent and which independent; mathematically, there's no easy way around that.

Hope this helps,

Erik Postma
Maplesoft.

Or, as a slight improvement, in this case you can also use ?LinearFit , since the expression is linear in the parameters (p[1] and p[2]). That will generally lead to a more accurate result and use less time. In fact, if you use ?Statistics[Fit] , then Maple will decide whether LinearFit is applicable or not; you could call it like this, for example:

independentDataMatrix := <<IC50> | <REDOX>>;
H := Statistics[Fit](p1 * x1 + p2 * x2 + p3, independentDataMatrix, LIPO, [x1, x2]);

You can add the output=solutionmodule option, like HerClau did above, if you like. Otherwise it will just give you the fitted function.

You need to make a decision which quantities you consider dependent and which independent; mathematically, there's no easy way around that.

Hope this helps,

Erik Postma
Maplesoft.

Hi emoh,

I think this is a case where we could possibly do better with documentation, explaining what the intent is, but I think that it does perform as intended. You wrote:

About VolumeOfRevolution it says:

"The VolumeOfRevolution(f(x), x=a..b) command returns the volume of revolution of the expression f(x) from a to b."

Which is not true.

If you plot
VolumeOfRevolution(2*x^2,x=0..2,axis=vertical,'output'='plot), the program instead returns the volume of revolution of 2*x^2 from y=0..8.

However, if you leave out the axis option (axis = vertical) or replace it with the default, axis = horizontal, then I think you get what you would like to get. (If not, feel free to reply again and correct me.) The axis option is explained further down on that page, ?Student,Calculus1,VolumeOfRevolution , as follows:

axis = horizontal or vertical Whether the expression is rotated horizontally or vertically. By default, the rotation is horizontal.

The explanation here is not very good - what is meant is that the axis of rotation is horizontal or vertical, respectively, not the direction in which the plot is swept around. (As an aside, it's also of course not the expression that is rotated, but the curve.) So if you add the axis=vertical option, then the curve which runs from x=0 to x=2 (and from y=0 to y=8) is rotated around the y-axis to obtain the "outer surface" of the volume described. Indeed, if you plot(2*x^2, x=0..2) and mentally rotate the resulting curve around either the horizontal axis (y=0) or the vertical axis (x=0), you get what I get if I call VolumeOfRevolution(2*x^2, x=0..2, axis=horizontal, output=plot) or the same with axis = vertical, respectively.

For the second point you make, I would say the same thing is happening: what you describe indicates to me that you want to rotate around the horizontal axis, not the vertical one. For the tutor, this may be what's going wrong again. Then you add the distancefromaxis=8 option; that again would make most sense to me in case axis = horizontal, because then the used axis of rotation would be the horizontal line y=8, which would touch the curve you are defining. Otherwise you would rotate around x=8, sweeping the curve in a big circle to (x, y) = (16, 0) and around.

Finally, for your example with computation, I get the answer Pi*128/3 when computing by hand for rotation around x=6: rotating around x=6 means that we need to view everything as a function of y, which means that there are two parts: from x=0 to x=6-y (for y=4..6) and from x=6-y to x=6 (for y=0..4). Since we're rotating around a shifted axis, we need to correct for that and obtain:

Pi * (int((6 - 0)^2 - (6 - (6 - y))^2, y=4..6) + int((6 - (6 - y))^2 - (6 - 6)^2, y = 0..4))

which simplifies to

Pi * (int(36 - y^2, y=4..6) + int(y^2, y=0..4))

and both by hand and using Maple that gives Pi * 128/3.

Your result by hand and the accompanying formula I could only obtain with different parameters: rotating x = 6-y, taken between y=0..4 (or equivalently, x=2..6) around x=6.

I hope this may clear things up a little bit,

Erik Postma
Maplesoft.

Hi emoh,

I think this is a case where we could possibly do better with documentation, explaining what the intent is, but I think that it does perform as intended. You wrote:

About VolumeOfRevolution it says:

"The VolumeOfRevolution(f(x), x=a..b) command returns the volume of revolution of the expression f(x) from a to b."

Which is not true.

If you plot
VolumeOfRevolution(2*x^2,x=0..2,axis=vertical,'output'='plot), the program instead returns the volume of revolution of 2*x^2 from y=0..8.

However, if you leave out the axis option (axis = vertical) or replace it with the default, axis = horizontal, then I think you get what you would like to get. (If not, feel free to reply again and correct me.) The axis option is explained further down on that page, ?Student,Calculus1,VolumeOfRevolution , as follows:

axis = horizontal or vertical Whether the expression is rotated horizontally or vertically. By default, the rotation is horizontal.

The explanation here is not very good - what is meant is that the axis of rotation is horizontal or vertical, respectively, not the direction in which the plot is swept around. (As an aside, it's also of course not the expression that is rotated, but the curve.) So if you add the axis=vertical option, then the curve which runs from x=0 to x=2 (and from y=0 to y=8) is rotated around the y-axis to obtain the "outer surface" of the volume described. Indeed, if you plot(2*x^2, x=0..2) and mentally rotate the resulting curve around either the horizontal axis (y=0) or the vertical axis (x=0), you get what I get if I call VolumeOfRevolution(2*x^2, x=0..2, axis=horizontal, output=plot) or the same with axis = vertical, respectively.

For the second point you make, I would say the same thing is happening: what you describe indicates to me that you want to rotate around the horizontal axis, not the vertical one. For the tutor, this may be what's going wrong again. Then you add the distancefromaxis=8 option; that again would make most sense to me in case axis = horizontal, because then the used axis of rotation would be the horizontal line y=8, which would touch the curve you are defining. Otherwise you would rotate around x=8, sweeping the curve in a big circle to (x, y) = (16, 0) and around.

Finally, for your example with computation, I get the answer Pi*128/3 when computing by hand for rotation around x=6: rotating around x=6 means that we need to view everything as a function of y, which means that there are two parts: from x=0 to x=6-y (for y=4..6) and from x=6-y to x=6 (for y=0..4). Since we're rotating around a shifted axis, we need to correct for that and obtain:

Pi * (int((6 - 0)^2 - (6 - (6 - y))^2, y=4..6) + int((6 - (6 - y))^2 - (6 - 6)^2, y = 0..4))

which simplifies to

Pi * (int(36 - y^2, y=4..6) + int(y^2, y=0..4))

and both by hand and using Maple that gives Pi * 128/3.

Your result by hand and the accompanying formula I could only obtain with different parameters: rotating x = 6-y, taken between y=0..4 (or equivalently, x=2..6) around x=6.

I hope this may clear things up a little bit,

Erik Postma
Maplesoft.

Ah, now I think I understand - I wasn't really focusing on the volume at all, just on the surface. And the surface that was generated was, I think, the correct one, but it is not really made clear what bounds the actual volume.

Of course one can do almost exactly what you did but "say it differently", somewhat more in line with the original poster's question, as follows:

with(Student[Calculus1]):
Q1 := VolumeOfRevolution(0, sqrt(x), x=0..4, axis=vertical, output=plot): Q1;
Q2 := SurfaceOfRevolution(4, y=0..2, axis=horizontal, output=plot, surfaceoptions=[color=blue]): Q2;
plots[display](Q1, plottools[rotate](Q2, 0, Pi/2, 0), title="All together now...");

Best,

Erik Postma
Maplesoft.

Ah, now I think I understand - I wasn't really focusing on the volume at all, just on the surface. And the surface that was generated was, I think, the correct one, but it is not really made clear what bounds the actual volume.

Of course one can do almost exactly what you did but "say it differently", somewhat more in line with the original poster's question, as follows:

with(Student[Calculus1]):
Q1 := VolumeOfRevolution(0, sqrt(x), x=0..4, axis=vertical, output=plot): Q1;
Q2 := SurfaceOfRevolution(4, y=0..2, axis=horizontal, output=plot, surfaceoptions=[color=blue]): Q2;
plots[display](Q1, plottools[rotate](Q2, 0, Pi/2, 0), title="All together now...");

Best,

Erik Postma
Maplesoft.

I just talked to one of our control theory experts. He says: "In fact, since his AR model is linear, the Kalman filter gives the best possible solution among all estimators as long as the random term is Gaussian.  (The Kalman filter is essentially applied to estimation problems with Gaussian noise). If the noise is not Gaussian the Kalman filter is still usable but it is suboptimal. "

If you would happen to have access to the Control Design Toolbox, you could use the ?ControlDesign,Kalman procedure. That would solve your potential implementation problems.

Hope this helps,

Erik Postma
Maplesoft.

I just talked to one of our control theory experts. He says: "In fact, since his AR model is linear, the Kalman filter gives the best possible solution among all estimators as long as the random term is Gaussian.  (The Kalman filter is essentially applied to estimation problems with Gaussian noise). If the noise is not Gaussian the Kalman filter is still usable but it is suboptimal. "

If you would happen to have access to the Control Design Toolbox, you could use the ?ControlDesign,Kalman procedure. That would solve your potential implementation problems.

Hope this helps,

Erik Postma
Maplesoft.

Right. I entered the coefficients that I saw in that figure into Maple just to see what happens. First of all, I think there are probably more significant digits than are shown in the images, because the fit was not very good. (The expression I ended up with was (0.841e-1+25.8927*x+3.5311*x^2-1.4071*x^3-.2573*x^4+0.42419e-3*x^5+0.15e-2*x^6+0.33739e-4*x^7)/(1+.1668*x-0.517e-1*x^2-0.117e-1*x^3-0.1981e-3*x^4+0.66445e-4*x^5+0.2534e-5*x^6).)

However, irrespective of these issues, the interpolant that MathCad finds has two asymptotes in the range that you're interested in: one at (about) 4.15 and one at (about) 13.0. You can see that in the following graph:

I would think that this is a more serious issue than the size of the residuals - you may have a good fit at the data points that you specify, but in between those data points the interpolant goes all the way to infinity and back !

So I would still suggest the use of an interpolant that has only one root - possibly just (x - a). In order to get a smaller residual, you can use some of the coefficients that were "freed up" by not using them in the denominator to increase the degree of the numerator. For example, you could use a 12th degree polynomial divided by (x - a). This gives a standard deviation of the residual of 0.12.

(Note that reducing the residual to these sizes is only useful if your data is actually that accurate - but I think it is indeed quite accurate, right?)

Hope this helps,

Erik Postma
Maplesoft.

Right. I entered the coefficients that I saw in that figure into Maple just to see what happens. First of all, I think there are probably more significant digits than are shown in the images, because the fit was not very good. (The expression I ended up with was (0.841e-1+25.8927*x+3.5311*x^2-1.4071*x^3-.2573*x^4+0.42419e-3*x^5+0.15e-2*x^6+0.33739e-4*x^7)/(1+.1668*x-0.517e-1*x^2-0.117e-1*x^3-0.1981e-3*x^4+0.66445e-4*x^5+0.2534e-5*x^6).)

However, irrespective of these issues, the interpolant that MathCad finds has two asymptotes in the range that you're interested in: one at (about) 4.15 and one at (about) 13.0. You can see that in the following graph:

I would think that this is a more serious issue than the size of the residuals - you may have a good fit at the data points that you specify, but in between those data points the interpolant goes all the way to infinity and back !

So I would still suggest the use of an interpolant that has only one root - possibly just (x - a). In order to get a smaller residual, you can use some of the coefficients that were "freed up" by not using them in the denominator to increase the degree of the numerator. For example, you could use a 12th degree polynomial divided by (x - a). This gives a standard deviation of the residual of 0.12.

(Note that reducing the residual to these sizes is only useful if your data is actually that accurate - but I think it is indeed quite accurate, right?)

Hope this helps,

Erik Postma
Maplesoft.

I'm not an expert on this, but the system that you posted looks like it might well be tackled with a tool from control engineers called a Kalman filter. That's a technique that also works online (that is, you get a useful estimation after a few data points that is improved as you continue) but it works just as well offline (where you just need the best possible estimate given all data points at once). Moreover, it is optimal in some way.

Only works if the "update function" is linear, though. (Well, there are similar approaches to nonlinear systems, but you lose the guarantees that it's optimal and even that it will converge at all, I believe.)

Hope this helps,

Erik Postma
Maplesoft.

PS: On an unrelated note, an optimization of the code you showed above. You only need to call Fit once:

w1, w2 := op(Fit(B2*x, XX, YY, x, output = [parametervector, standarderrors])):

will do the trick.