What is the best approach to calculating function limits?

Question

I am planning to develop an application that can calculate the limit of an expression (function) given by the user.

I already have a functional expression evaluator, which will definitely come in handy. My idea is to calculate it like this: I give the parameter a few values that get closer and closer to the point, but doesn't reach the point. At the end, I see if the difference between two consecutive results gets smaller or bigger closer or further from 0. If it gets smaller closer, this means that the limit is finite, infinite otherwise. After that, it's easy to approximate a result.

Is there a better approach to this? Or should I go with this one?

The application will accept functions that contain these mathematical operators: +,-,*,/,%,^, functions (like floor, logarithms, trigonometry), and 3 condition functions (abs, min, max). So, for example a function which has a specific value for integer values, and another value for non-integer values is not accepted.

Answer 1

Although it is hard to find a limit directly, you can estimate them pretty well in C# using the double.Epsilon constant to find an extremely close value. This field is a constant for the smallest number that can be represented by a double floating point number, so we'll use it to find an extremely close value to the x value and plug that into a function instead.

Here is an example:

using System;

namespace LimitsTest
{
    class Functions
    {
        public static double OneDividedByX(double x)
        {
            return 1 / x;
        }

        public static double XDividedByX(double x)
        {
            return x / x;
        }

        public static double OneDividedByXSquared(double x)
        {
            return 1 / Math.Pow(x, 2);
        }

        public static double XDividedByXPlusTwo(double x)
        {
            return x / (x + 2);
        }
    }

    class Program
    {
        static double LeftLimit(double x, Func<double, double> function)
        {
            return function(x - double.Epsilon);
        }

        static double RightLimit(double x, Func<double, double> function)
        {
            return function(x + double.Epsilon);
        }

        static double Limit(double x, Func<double, double> function)
        {
            double right = LeftLimit(x, function);
            double left = RightLimit(x, function);

            return (right == left) ? right : double.NaN;
        }

        static void ShowLimits(double x, string description, Func<double, double> function)
        {
            Console.WriteLine("{0} at {1}: {2}", description, x, function(x));
            Console.WriteLine("Limit of {0} as approached from the left of {1}: {2}",
                description, x, LeftLimit(x, function));
            Console.WriteLine("Limit {0} as approached from the right of {1}: {2}",
                description, x, RightLimit(x, function));
            Console.WriteLine("Limit of {0} as approached from {1}: {2}",
                description, x, Limit(x, function));
            Console.WriteLine();
        }

        static void Main(string[] args)
        {
            ShowLimits(0, "1 / x", Functions.OneDividedByX);
            ShowLimits(0, "x / x", Functions.XDividedByX);
            ShowLimits(0, "1 / x^2", Functions.OneDividedByXSquared);
            ShowLimits(-2, "x / (x + 2)", Functions.XDividedByXPlusTwo);

            Console.ReadLine();
        }
    }
}

Answer 2

The way I did this is first apply Taylor to find an estimate of the coefficients for the function given but shifts the function to left by the limit x. Then from their set coefficient that are nan to 0. Then just calculate the function at the x value given.

#include <iostream>
#include <vector>
#include <math.h>

double getDerivitive(unsigned derivitive, double x, double cx, const std::function<double(double)> &f) {
    if (derivitive == 0) return f(x);

    double newCx = derivitive * cx;

    return (getDerivitive(derivitive - 1, x+newCx, cx, f) - getDerivitive(derivitive - 1, x-newCx, cx, f)) / (2*newCx);
}
std::vector<double> getCoefficient(unsigned precision, const std::function<double(double)> &f) {

    std::vector<double> coeffs(precision, 0.0);
    double iFactorial = 1;
    for (unsigned i = 0; i < precision; i++) {

        coeffs[i] = getDerivitive(i, 0.0, 0.01, f) / iFactorial;

        iFactorial *= (i+1);
    }

    return coeffs;
}

double getValueAt(double x, const std::vector<double> &coeffs) {

    double total = 0;

    double xToTheI = 1;
    for (unsigned i = 0; i < coeffs.size(); i++) {

        total += coeffs[i] * xToTheI;

        xToTheI *= x;
    }

    return total;
}

double getLimit(double x, const std::function<double(double)> &f) {

    std::function<double(double)> newFunc = [&](double atX) {
        return f(x + atX);
    };

    std::vector<double> coefficients = getCoefficient(10, newFunc);

    if (isnan(coefficients[0])) {
        coefficients[0] = newFunc(0.00001);

        if (abs(coefficients[0]) < 0.0001) {
            coefficients[0] = 0;
        }
    }

    for (unsigned i = 0; i < coefficients.size(); i++) {
        if (isnan(coefficients[i])) {
            coefficients[i] = 0;
        }
    }

    return getValueAt(0, coefficients);
}

int main(int argc, const char * argv[]) {


    std::function<double(double)> cosFTest = [](double x) {
        return cos(x);
    };

    std::cout << "Limit of cos(0)  ==  " << getLimit(0, cosFTest) << std::endl;

    std::function<double(double)> sinOverXTest = [](double x) {
        return sin(x)/x;
    };

    std::cout << "Limit of sin(x)/x  ==  " << getLimit(0, sinOverXTest) << std::endl;

    std::function<double(double)> polynomial = [](double x) {
        return (x*x) / x;
    };

    std::cout << "Limit of x^2 / x  ==  " << getLimit(0, polynomial) << std::endl;

    std::function<double(double)> funcTestInfinity = [](double x) {
        return 1.0/x;
    };

    std::cout << "Limit at 1/x apraoches infinity ==  " << getLimit(INFINITY, funcTestInfinity) << std::endl;

    return 0;
}

Answer 3

mathematically you could use the Differential calculus. Like that you just have to implement the differential rules and you dont have to use brute force

but the comment you got with wolframalpha is great: example: thats exactly what you need

Answer 4

This answer is more mathematics than programming, but it shows you why you can't do what you are looking for (unless you add more information about the functions).

We define f as follows:

• f(x) = 0 if x is a rational number
• f(x) = 1 if x is not a rational number

This function doesn't even have a limit at any point, but if you use the method you specified, then for what ever float (or double) number you use, f(x) will be 0, so the difference will be 0.

Your method will say that the limit of f at the point 5 for example is 0, but f doesn't have a limit at all.

Answer 5

I was just going to write a brief comment, but the more I thought about it, the more I came to the conclusion that the answer to your question is that you should not go with your suggestion. There is a circularity in your proposal: to find the limit of a function you propose to provide your system with a series of inputs whose evaluation tends towards the (supposed) limit of the function. How would you start without having an idea of what the limit is ? Or do you want to implement a program which finds the limits of functions that you already know ? A sort of artificial idiot who will always know less than you ?

As to what you should do, first realise that this is actually quite a tricky function to implement and that there will be no general solution that works for all functions, not even for all the nice, well-behaved functions. You may find it instructive to see how Sage (http://www.sagemath.org/), or any other open-source computer algebra system you know about, does this.