Home > Mobile >  How do I implement the numerical differentiation (f'(x) = f(x h)-f(x)/ h
How do I implement the numerical differentiation (f'(x) = f(x h)-f(x)/ h

Time:04-23

2nd task:

For a function f : R^n → R the gradient at a point ~x ∈ R^n is to be calculated: - Implement a function CMyVector gradient(CMyVector x, double (*function)(CMyVector x)), which is given in the first parameter the location ~x and in the second parameter the function f as function pointer in the second parameter, and which calculates the gradient ~g = grad f(~x) numerically by gi = f(x1, . . . , xi-1, xi h, xi 1 . . . , xn) - f(x1, . . . , xn)/h to fixed h = 10^-8.

My currently written program:

Header
#pragma once
#include <vector>
#include <math.h>
class CMyVektor
 {
   private:
        /* data */
       int Dimension = 0;
        std::vector<double>Vector;
    public:
        CMyVektor();
        ~CMyVektor();
        //Public Method
        void set_Dimension(int Dimension /* Aktuelle Dim*/);
        void set_specified_Value(int index, int Value);
        double get_specified_Value(int key);
        int get_Vector_Dimension();
        int get_length_Vektor();
        double& operator [](int index);
        string umwandlung()
 };
 CMyVektor::CMyVektor(/* args */)
 {
     Vector.resize(0, 0);
 }
 CMyVektor::~CMyVektor()
 {
     for (size_t i = 0; i < Vector.size(); i  )
     {
         delete Vector[i];
     }
 }
void CMyVektor::set_Dimension(int Dimension /* Aktuelle Dim*/)
{
        Vector.resize(Dimension);
};
void CMyVektor::set_specified_Value(int index, int Value)
{
    if (Vector.empty()) 
    {
        Vector.push_back(Value);
    }
    else {
        Vector[index] = Value;
    }
};
double CMyVektor::get_specified_Value(int key)
{
    // vom intervall anfang - ende des Vectors
    for (unsigned i = 0; i < Vector.size(); i  )
    {
        if (Vector[i] == key) {
            return Vector[i];
        }
    }
};
int CMyVektor::get_Vector_Dimension()
{
    return Vector.size();
};
// Berechnet den Betrag "länge" eines Vectors.
int CMyVektor::get_length_Vektor() 
{
    int length = 0;
    for (size_t i = 0; i < Vector.size(); i  )
    {
        length  = Vector[i]^2
    }
    return sqrt(length);
}

// [] Operator überladen
double& CMyVektor::operator [](int index)
{
    return Vector[index];
}
main.cpp

#include <iostream>
#include "ClassVektor.h"

using namespace std;

CMyVektor operator (CMyVektor a, CMyVektor b);
CMyVektor operator*(double lambda, CMyVektor a);
CMyVektor gradient(CMyVektor x, double (*funktion)(CMyVektor x));

int main() {

    CMyVektor V1;
    CMyVektor V2;


    CMyVektor C;
    C.set_Dimension(V1.get_length_Vector());
    C= V1   V2;
    std::cout << "Addition : "<< "(";;
    for (int i = 0; i < C.get_length_Vector(); i  )
    {
            std::cout  << C[i] << " ";
    }
    std::cout << ")" << endl;
    C = lamda * C;
    std::cout << "Skalarprodukt: "<< C[0]<< " ";

}
// Vector Addition
CMyVektor operator (CMyVektor a, CMyVektor b)
{
    int ai = 0, bi = 0;
    int counter = 0;
    CMyVektor c;
    c.set_Dimension(a.get_length_Vector());
    // Wenn Dimension Gleich dann addition
    if (a.get_length_Vector() == b.get_length_Vector())
    {
        while (counter < a.get_length_Vector())
        {
            c[counter] = a[ai]   b[bi];
            counter  ;
        }
        return c;
    }

}
//Berechnet das Skalarprodukt
CMyVektor operator*(double lambda, CMyVektor a) 
{
    CMyVektor c;
    c.set_Dimension(1);
    for (unsigned i = 0; i < a.get_length_Vector(); i  )
    {
        c[0]  = lambda * a[i];
    }
    return c;
}

/*
* Differenzenquotient : (F(x0 h) F'(x0)) / h
* Erster Parameter die Stelle X - Zweiter Parameter die Funktion
* Bestimmt numerisch den Gradienten.
*/
CMyVektor gradient(CMyVektor x, double (*funktion)(CMyVektor x))
{

}

My problem now is that I don't quite know how to deal with the CMyVector gradient(CMyVector x, double (*function)(CMyVector x)) function and how to define a function that corresponds to it.

I hope that it is enough information. Many thanks.

CodePudding user response:

The function parameter is the f in the difference formula. It takes a CMyVector parameter x and returns a double value. You need to supply a function parameter name. I'll assume func for now.

I don't see a parameter for h. Are you going to pass a single small value into the gradient function or assume a constant?

The parameter x is a vector. Will you add a constant h to each element?

This function specification is a mess.

Function returns a double. How do you plan to turn that into a vector?

No wonder you're confused. I am.

Are you trying to do something like this?

CodePudding user response:

You are given a function signature

CMyVector gradient(CMyVector x, double (*function)(CMyVector x))

Without knowing the exact definition I will assume, that at least the basic numerical vector operations are defined. That means, that the following statements compile:

CMyVector x {2.,5.,7.};
CMyVector y {1.,7.,4.};
CMyVector z {0.,0.,0.};
double a = 0.;
// vector addition and assigment
z = x   y;
// vector scalar multiplication and division
z = z * a;
z = x / 0.1;

Also we need to know the dimension of the CMyVector class. I assumed and will continue to do so that it is three dimensional.

The next step is to understand the function signature. You get two parameters. The first one denotes the point, at which you are supposed to calculate the gradient. The second is a pointer to the function f in your formula. You do not know it, but can call it on a vector from within your gradient function definition. That means, inside of the definition you can do something like

double f_at_x = function(x);

and the f_at_x will hold the value f(x) after that operation.

Armed with this, we can try to implement the formula, that you mentioned in the question title:

CMyVector gradient(CMyVector x, double (*function)(CMyVector x)) {
    double h = 0.001;
    
    // calculate first element of the gradient
    CMyVector e1 {1.0, 0.0, 0.0};
    double result1 = ( function(x   e1*h) - function(x) )/h;

    // calculate second element of the gradient
    CMyVector e2 {0.0, 1.0, 0.0};
    double result2 = ( function(x   e2*h) - function(x) )/h;
    
    // calculate third element of the gradient
    CMyVector e3 {0.0, 0.0, 1.0};
    double result3 = ( function(x   e3*h) - function(x) )/h;

    // return the result
    return CMyVector {result1, result2, result3};
}

There are several thing worth to mention in this code. First and most important I have chosen h = 0.001. This may like a very arbitrary choice, but the choice of the step size will very much impact the precision of your result. You can find a whole lot of discussion about that topic here. I took the same value that according to that wikipedia page a lot of handheld calculators use internally. That might not be the best choice for the floating point precision of your processor, but should be a fair one to start with.

Secondly the code looks very ugly for an advanced programmer. We are doing almost the same thing for each of the three dimensions. Ususally you would like to do that in a for loop. The exact way of how this is done depends on how the CMyVector type is defined.

CodePudding user response:

Since the CMyVektor is just rewritting the valarray container, I will directly use the valarray:

#include <iostream>
#include <valarray>

using namespace std;

using CMyVektor = valarray<double>;
CMyVektor gradient(CMyVektor x, double (*funktion)(CMyVektor x));

const double h = 0.00000001;

int main()
{
    // sum(x_i^2   x_i)--> gradient: 2*x_i   1
    auto fun = [](CMyVektor x) {
        double res;
        for(auto i: x) res  = i*i   i;
        return res ;};
     
    CMyVektor d = gradient(CMyVektor{1,2,3,4,5}, fun);
    for (auto i: d) cout << i<<' ';
    return 0;
}

CMyVektor gradient(CMyVektor x, double (*funktion)(CMyVektor x)){
    CMyVektor grads(x.size());
    CMyVektor pos(x.size());
    for (int i = 0; i<x.size(); i  ){
        pos[i] = 1;
        grads[i] = (funktion(x   h * pos) - funktion(x))/ h;
        pos[i] = 0;
    }
    return grads;
}

The prints out 3 5 7 9 11 which is what is expected from the given function and the given location

  • Related