Find the parameter $a$ and the coordinate $x$ of a function so that the function's maximum has a given value

Question

I have implemented in C++ the code for a quasi-Newton method given by

$x_{n+1} = x_n - \dfrac{\delta f(x_n)}{f(x_n+\delta) - f(x_n)}$

for a computational project in which $f$ is a function given by:

$f(x) = 2e^{-ax}(a\sin(3x) + \cos(3x))$

One of the exercises asks to find the parameter $a$ and the coordinate $x$ so that the maximum value of the function $f$ is 3. In the problem is left as a suggestion that we use the above quasi-Newton method implemented to find the value of $x$ and the bisection method to find $a$. Any tip on how to do this?

Thanks,

Patrick

Answer 1

Your Quasi-Newton method is not correct, you need a $-$ after $x_n$. Note also that if $x$ is a local maximum, then $f'(x)=0$ and to solve this with your Quasi-Newton you don't have to provide $f''(x)$ analytically, as it's already estimated "in place" with finite differences ( i.e. $f'(x_0)= \frac{f(x_0+\delta) - f(x_0)}{\delta} + O(\delta)$).

After a quick test, I obtained $x = -0.09679574815$, $a= -6.209538561$. A plot confirms that the local maximum is $3$ at point $x$.

Another found solution is $x=-2.052921249$ and $a=0.1898554466$.

I applied "externally" the bisection method for the parameter $a$ and inside of it I used the quasi-Newton method in order to solve for $x$ the equation $f_a'(x)=0$ for different parameters $a$. As $(x,a)$ are coupled, you can't split the two methods.

I wrote the following C++ snippet: here a1 and a2 are the leftmost and rightmost point of the starting interval in the bisection routine.

    #include <cmath>
    #include <iostream>
    #include <iomanip>      // std::setprecision
    
    constexpr double tol = 1e-12;
    constexpr double delta = 1e-5;
    double f(double x, double a){
        return 2*exp(-a*x)*(a*sin(3*x)+cos(3*x));
    }
    
    double fprime(double x, double a){
        return exp(-a*x)*(4*a*cos(3*x) - 2*(a*a+3)*sin(3*x));
    }
    
    double quasi_newton(double f(double ,double),double x0,double a){
        double res =  -(delta*f(x0,a) )/(f(x0+delta,a)-f(x0,a));
        while (std::fabs(res)>tol) {
            x0+=res;
            res = -(delta*f(x0,a))/(f(x0+delta,a)-f(x0,a));
        }
        return res+x0;
    }
    
    
    
    int main(){
        double x0=0.5;
        int it{0};
        double a1{-10.0};
        double a2{10.0};
        double c{0};
        
        double xc,gc,ga1,ga2,xa1,xa2;
        while (std::fabs(a1-a2)>tol) {
            std::cout<< std::fabs(a1-a2) <<std::endl;
            c =0.5*(a1+a2);
            xc = quasi_newton(fprime,x0,c);
            xa1 = quasi_newton(fprime,x0,a1);
            xa2 = quasi_newton(fprime,x0,a2);
            
            gc = f(xc,c)-3.0;
            ga1 = f(xa1,a1)-3.0;
            ga2 = f(xa2,a2)-3.0;
    
            if (gc*ga1<0) {
                a2 = c;
            }
            if (gc*ga2<0) {
                a1=c;
            }
            
        }
        std::cout <<"Point: " << std::setprecision(10) << xc << ", Value of a: " <<std::setprecision(10)<< c << ", Check:" << std::setprecision(10) << f(xc,c) << std::endl;
        return 0;
        
    }