Testing AI Tools with a Integral Calculus Task

All this AI buzz got me curious, so I decided to test some of the products out there by using them to solve for an antiderivative.

This problem looks at one aspect of my previous work involving elliptical motion. See the diagram below.

When a body's motion is caused by a centrally-directed force, for example, one body revolving around another, it sweeps out a constant area per unit time. In the case of an ellipse, this motion can be described in polar coordinates with origin based on the major axis at the focus of the motion (e, 0). The radial vector that sweeps out the motion, r_η, is a function of the angle eta (η) measured relative to the major axis:

To simplify work, a will be set equal to 1 (it can be scaled appropriately as needed; the hard part is solving for position as a function of time). In addition, the same values will be used as in previous posts: T = 2*π time units and e = 0.786151377748.

The area of an ellipse, A_ellipse, is π a b.
In the present case, since a = 1, this area is

A_ellipse = π √(1 - e²)

One quarter of this area is

A_{quarter-ellipse} = π √(1 - e²) / 4

See the following diagram.

In this post, I am going to find the area when r_η coincides with the semi-latus rectum; in other words, when η = π/2 and r_η = 1 - e². See the following diagram.

Note that this area is not the same as the area for one-quarter of the ellipse--unless e = 0 (i.e., it's a circle).

This area can be calculated in at least four ways. Ideally, I'd like to find an expression that does not involve a trigonometric function.

1) Kepler's original process of geometry. Consider the area swept out by the superposing radius, remove the area of the triangle eOP, then scale that area by √(1 - e²):

A_tp = √(1 - e²) [0.5 * a² * E - 0.5 * base * height]
    = 0.5 * √(1 - e²) [1 * arccos(e) - e * √(1 - e²)]
    = 0.5 * √(1 - e²) [arccos(e) - e * √(1 - e²)]
    = 0.0557377539689 units of area.

2) Knowing that equal area is swept out in equal time, dA/dt = constant:

dA/dt = A_ellipse / T = π √(1 - e²)/(2*π) = √(1 - e²)/ 2

The time at which r_η coincides with the semi-latus rectum is

t_p = E - e sin(E) = arccos(e) - e √(1 - e²) = π/2 - arcsin(e) - e √(1 - e²)

A_tp = dA/dt * t_p
    = 0.5 * √(1 - e²) * [π/2 - arcsin(e) - e √(1 - e²)]
    = 0.309016994381 * 0.180371160753
    = 0.0557377539689 units of area.

(Even though the expression with arccos is shorter, I prefer the one with arcsin because it explicitly shows how this quantity deviates from π/2.)

3) When η = π/2, the r_η vector is a vertical straight line, and standard integration with cartesian coordinates can be used.

The antiderivative of this function is

Evaluating this antiderivative at the endpoints to get the definite integral yields

Evaluated, A_tp = 0.0557377539689 units of area.
This number is the same as found above.
Note that this expression highlights the fact that
as e goes to 0, A_tp goes to A_{quarter-ellipse}.

4) A_tp can also be found by computing the definite integral of the original equation in polar coordinates. The area swept out by a vector defined in terms of a polar equation is given as follows:

In the present case, the definite integral is

The 0.5 and (1 - e²)² factors are constant, so they can be taken out of the integrand, leaving an integral of the form

This integral is not straightforward to solve. Its antiderivative can be derived, but I looked it up in a table of antiderivatives (the good old CRC Handbook):

The antiderivative of the integral in this antiderivative, call it I₂, can be expressed in two ways.

However, one of those expressions includes the term √(b² - a²). But because in the present context (for elliptical motion 0 < e < 1) b is always less than a that would yield an imaginary result. So we are left with the following result:

Subbing in a = 1 and b = e the definite integral becomes

Now the two constant factors that were removed, 0.5 and (1 - e²)², can be put back and the area of the region of interest is

Unfortunately, the expressions for A_tp seem to be getting more complicated--and there is still the presence of a trigonometric function.

Plugging in the numbers, A_tp = 0.0557377539689 units of area.
This number is the same as found above, so this number is verified.

Actually, a factor of √(1 - e²) can be taken out of the denominator to yield an expression very similar to those produced by the first three methods.

Here is where AI came in. I asked several of the AI programs two questions:

(i) what is the antiderivative of the original polar equation for area, and
(ii) what is the definite integral of the original polar equation for area.

Most of them did not provide an answer. Those that did were wrong.

Here is the response from grok:

And here is grok's response to my request to go directly to the definite integral:

Neither of these answers is confirmed by the work above. Looking at grok's first response, if the antiderivative is evaluated at the endpoints, it does not yield the same numeric values found above.

And grok's second response is wrong too.
Multiplying this result by (1 - e²)²/2 yields π √(1 - e²)/4,
which is A_{quarter-ellipse}--the incorrect area.

Perplexity did not get it.
ChatGPT did not get it.

Asking Google from its main search page did not yield a response--period. (For a couple weeks.)

However, asking directly on the Gemini home page did yield results--and they are correct:

This result can be simplified further:

Returning the factors 0.5 and (1 - e²)², the final expression for A_tp can be expressed as

This expression could be simplified further by factoring the (1 - e²) term out of the denominator, which would yield an expression identical to the ones above, but leaving the (1 - e)^3/2 term in the denominator does serve to remind us of one of Kepler's laws of celestial motion: the square of a planet's orbital period is proportional to the cube of the semi-major axis of its orbit. In other words, T ∝ a^3/2.

Working through the numbers again,

A_tp = 0.764064498063 * 0.072949016881
    = 0.0557377539674

So the results above have been further confirmed.





Copilot got it. Here's its result for the indefinite integral:

And here's its result for the definite integral:

Both these answers confirm previous work.


This post does not cover new ground, I haven't discovered anything new. But I have found some weaknesses of the current crop of AI products.

Kepler's Equation: A Listing of Some Derivatives

This post serves as a compilation of some derivatives for Kepler's Equation of Elliptical Motion. Nothing fancy here. Just posting some basic characteristics of the motion that might be useful later.

For a start, here are the M and E curves plotted on the same graph:

From Kepler's Equation of Elliptical Motion, the derivative of E with respect to (w.r.t.) t is found:

For the model I've been using throughout this work, ω = 1, so this derivative can be simplified further. In fact, I'm going to examine derivatives w.r.t. M instead of time. RETURNING TO A GENERIC FORM IN WHICH ω IS NOT 1, WE HAVE TO REMEMBER TO PUT ω BACK IN.

Here is the derivative of E w.r.t. M:

Here is a graph of dE/dM:

Here is the second derivative of E w.r.t. M:

Here is the third derivative of E w.r.t. M:

Each term in the numerator of this expression has a factor of e; when e is 0, the third derivative is 0.

This derivative also has a zero at a non-intuitive point. To find this root, rearrange this equation so that it is solely in terms of cos(E) or sin(E). The simpler approach is to express it in terms of cos(E):

In this form, it is a quadratic expression in terms of cos(E); solving it yields two roots. One root is complex, which is meaningless in the present context. The other root follows:

Numerically, this corresponds to:
cos(E) = 0.947351608247
E = 0.325935469297
sin(E) = 0.320195144173
M = 0.074213615557

I haven't seen these numbers up to now. (I wonder what their significance is in terms of the motion?)

Carrying on, here is the fourth derivative of E w.r.t. M:

The numerator of this expression is already in the form of a quadratic equation in terms of cos(E). Again, one root is complex, which is meaningless in the present context. The other root follows:

Numerically, this corresponds to:
cos(E) = 0.85192985408
E = 0.551136659451
sin(E) = 0.523655921123
M = 0.139463835594

These derivatives may be useful in creating a Taylor polynomial and/or creating interpolating functions over certain intervals.

One comment I would like to make immediately is that, in using the quotient rule to perform the differentiation, the term (1 - e cos(E)) always appears in the numerator before simplification--it's not present in just the denominator. Simplifying the numerator may reveal some insight into the motion; however, if the derivatives are going to be used in a Taylor polynomial, and computed numerically, it may be a better idea to leave the (1 - e cos(E)) term as is. After all, this term has to be computed for the denominator anyhow. The better course of action may be to compute the term (1 - e cos(E)) once, save its value in a variable, and re-use it as needed. Just a thought. It appears in all the derivatives that follow too.

Some features to note:
(i) When the motion crosses the y-axis (i.e., E = π/2 rad),
dE/dt = ω (i.e., dE/dt is parallel to M; i.e., dE/dM = 1);
(ii) dE/dM is a maximum at t = 0, at which time it takes the value 1/(1 - e);
(iii) dE/dM is a minimum when M = E = π, at which time it takes the value 1/(1 + e).
(iv) from the second derivative of E w.r.t. M, the point of inflection for the E curve occurs when sin(E) is 0.

Carrying on, here are the cos(M) and cos(E) curves plotted on the same graph:

The red line is cos(E) and the blue line is cos(M).

The derivative of cos(E) w.r.t. M is found as

Here are plots of the first derivatives of the cosine curves on the same graph:

The second derivative of cos(E) w.r.t. M is

Here are plots of the second derivatives of the cosine curves on the same graph:

Note that the second derivative of cos(E) takes a zero value
when cos(E) = e; in other words, the point of inflection occurs when the motion is at the semi-latus rectum (M = 0.180371160753 rad).
Also note that when elliptical motion reduces to circular motion,
(i.e., e = 0), the point of inflection occurs at the usual times:
every time cos(M) = 0 (i.e., M = π /2, 3 π /2, etc.).

Now repeat the process for the sine curves. Here are the sin(M) and sin(E) curves plotted on the same graph:

The red line is sin(E) and the blue line is sin(M).

The derivative of sin(E) w.r.t. M is

Here are plots of the first derivatives of the sine curves on the same graph:

The second derivative of sin(E) w.r.t. M is

Note that this derivative is similar to the second derivative of E w.r.t. M; it's just missing the scaling factor e.

Here are plots of the second derivatives of the sine curves on the same graph:

Note that the second derivative of sin(E) takes a zero value
when sin(E) = 0; in other words, points of inflection occur each time M = 0, π, 2π, etc.
Again note that when elliptical motion reduces to circular motion,
(i.e., e = 0), the point of inflection occurs at the usual times:
every time sin(M) = 0 (i.e., M = 0, π , 2 π , etc.).

Kepler's Equation: Eccentric Anomaly Values at the Quarter-period - Part 2

I am circling back around to post some additional information about elliptical motion at the quarter period.

According to Kepler's Equation for Elliptical Motion,

E = M + e sin(E)                     Equation 1

At the quarter period, the Mean Anomaly, M, becomes π/2 and this equation becomes

E = π/2 + e sin(E)                     Equation 2

A plot of sin(E) at T/4 points for e going from 0 to 1 follows.

This graph, and the data that produced it, are posted on the Desmos website: Plot of sin(E) at the Quarter Period

These data points were computed numerically.
The first column, labelled x₁, contains the e values.
The second column is the Eccentric Anomaly, E, at T/4.
The third column is sin(E) at T/4.
This data is available for use, for anybody interested. Feel free to edit it, or plot functions using the data. However, if you want to save your work, you'll have to create your own account first (it's free.)

Equation 2 can be reduced to

sin(E) = cos[e sin(E)]                     Equation 3

or

x = cos(e x)                     Equation 4

where x = sin(E).

Note that when e = 1, Equation 4 reduces to the familiar form of the definition of the Dottie Number: x = cos(x). In other words, the Dottie Number is just one point on the curve that solves Equation 4.

Now consider the first derivative of E with respect to (w.r.t.) e. Equation 1 is differentiated w.r.t. e to yield

Equation 5 is used to differentiate sin(E) w.r.t. e:

or

A plot of d[sin(E)] / de at the quarter period follows:

This graph, and the data that produced it, are posted on the Desmos website: Plot of d[sin(E)]/de at the Quarter Period

Next, consider the second derivative of sin(E) w.r.t. e:

A plot of d²[sin(E)] / de² at the quarter period follows:

This graph, and the data that produced it, are posted on the Desmos website: Plot of d²[sin(E)]/de² at the Quarter Period

This plot takes a zero value at approximately x = 0.643638135454661. This means the plot for sin(E) values at the quarter period has a point of inflection. For 0 < e < 0.643638135454661 the plot for sin(E) is concave down. For 0.643638135454661 < e < 1 the plot for sin(E) is concave up.

At this point, d[sin(E)] / de = -0.33323545 and
d[cos(E)] / de = -0.5446936

Going through similar steps for cos(E), analogous to Equation 3, values for E_T/4 can be found as the solution to the equation

cos(E) = -sin[e sin(E)]

When e = 1, cos(E_T/4) = - sin(D), where D is the Dottie Number.

Here is a plot of cos(E) at T/4 points for e going from 0 to 1.

The first and second derivatives are

Plots of these derivatives are posted at the following URLs:

First Derivative of cos(E) w.r.t. e at the Quarter Period

Second Derivative of cos(E) w.r.t. e at the Quarter Period

Since the second derivative does not take a zero value over the range of interest ([0, 1]), the curve for cos(E) does not have a point of inflection over this range. Because the concavity of this curve does not change, perhaps it would be more effective to examine features of the cosine curve instead of the sine curve (e.g., a straight line bounding function would not have the risk of the function diverging away from it.)

Here is a plot of cos(E_T/4) with two simple bounding functions:

cos(E) at the Quarter Period with Two Simple Bounding Functions

Investigating the cos(E) curve, recall the equation stated above:

cos(E) = -sin[e sin(E)]

Also, consider the infinite series definition of the sine function:

In the present case, x = e sin(E). Notice that a factor of e is present in all terms.

So, cos(E) can be stated in the form

cos(E) = -e F(e)

where F(e) is some function to be determined.

The corresponding form of sin(E) is

sin(E) = √[1 - e² (F(e))²]

What does a plot of F(e) look like? Here is a plot of F(e) = -cos(E)/e:

Graph of -cos(E) divided by e at the Quarter Period

Note that F(e) goes to 1 as e goes to zero,
and it goes to sin(D) as e goes to 1.

In addition, the series representation indicates that the expression at the sin(E) = sin(M) point could be expressed in a similar form:

cos(E_sE=sM) = -sin[(e/2) sin(E)] = - (e/2) F(e/2)

These proposed forms fit the graphs seen so far.

Looking at other properties of E at the quarter period ...