The Newton Raphson method for p(x) = 100 .* sqrt(x) .* 0.99.^(sqrt(x)) - 0.5 .* x, problem with derivative and reasonable x value

Hi!

I'm very new to coding an MATLAB and only understand about half of what I'm doing and would really apricate some help.

I'm trying to solve

p (x) = 100 .* sqrt(x) .* 0.99.^(sqrt(x)) - 0.5 .* x

numerically. I have graphed it and found that the critical point x should be around 1600, se picture

Using the code:

x=[1598:1602];
p= 100 .* sqrt(x) .* 0.99.^(sqrt(x)) - 0.5 .* x;

figure
plot (x, p)
xlabel('Material (x)')
ylabel('Profit (p)')

Now I'm trying to find the exact x and I have found a lot of diffrent codes online to build the The Newton Raphson method, but are very confused as there is so many ways of doing it. I have tried a couple of diffrent code, but non have worked. Or they have given me values but the valus for x are nowhere near 1600. And I don't understand what is wrong, I'm thinking it might not be the code but the derivative I get wrong, or somthing.

This is one of the scripts I have tried for Newton method

% Newton's Method

f = @(x) 100 .* sqrt(x) .* 0.99.^(sqrt(x)) - 0.5 .* x;
fdiff = @(x) 50*(99/100)^(x^(1/2))*log(99/100) + (50*(99/100)^(x^(1/2)))/x^(1/2) - 1/2;

x0 = 1601;

iterations = 10;

xnext = x0;

for i=1:iterations
    xnext = xnext - f(xnext)/fdiff(xnext);
end

display("With " + iterations + " iterations of Newton's Method");
display(xnext);
display(abs(sqrt(2)-xnext));

And i have gotten the fdiff by doing the following in the command window

syms x
f=inline('100 .* sqrt(x) .* 0.99.^(sqrt(x)) - 0.5 .* x', 'x')

f =

     Inline function:
     f(x) = 100 .* sqrt(x) .* 0.99.^(sqrt(x)) - 0.5 .* x

diff(f(x),x)
 
ans =
 
50*(99/100)^(x^(1/2))*log(99/100) + (50*(99/100)^(x^(1/2)))/x^(1/2) - 1/2
 

If I try to run this I get

newtons_method

"With 10 iterations of Newton's Method"

xnext =

3.3402e-05 - 5.7903e-24i

1.4142

It's saying x is 3.3402e-05 - 5.7903e-24i

which isn't reasonable.

If I try this code wich is a code I found online and that I asked ChatGPT to help me taylor to my function (I know its not the most reliable, but at this point I don't know what to do).

function [root, number_of_iteration] = newton_profit(initial, tolerance, maxiteration)
    % Define the function
    f = @(x) 100 .* sqrt(x) .* 0.99.^(sqrt(x)) - 0.5 .* x;
    
    % Define the derivative of the function
    f_prime = @(x) 50 .* (0.99.^(sqrt(x))) .* (log(0.99) ./ (2 * sqrt(x)) + 1 ./ (2 * sqrt(x))) - 0.5;
    
    % Initial guess
    x = initial;
    number_of_iteration = 0;

    while abs(f(x)) > tolerance
        number_of_iteration = number_of_iteration + 1;
        
        if number_of_iteration > maxiteration
            error('Failed to converge, the maximum number of iterations is reached');
        end
        
        % Calculate the value of the function and its derivative
        sum = f(x);
        dif = f_prime(x);
        
        if dif == 0
            error('The derivative of the function is zero, pick another initial point and run again');
        end
        
        % Newton's formula
        x = x - sum / dif;
    end
    
    root = x;
end

% Example usage
initial_guess = 1600; % Start around your known value
tolerance = 1e-5; % Define the desired tolerance
maxiteration = 100; % Maximum number of iterations allowed

[root, num_iterations] = newton_profit(initial_guess, tolerance, maxiteration);
disp(['Root: ', num2str(root)]);
disp(['Number of Iterations: ', num2str(num_iterations)]);

I then get

Root: 7235.8192

Number of Iterations: 15

Witch also isn't reasonable.

Remeber the graph, it should def be around 1600

This is my first cours in math at university so I might also have missed something important in the math. But the answer should be around 1600 right? Can somone please help me to understand what is going wrong and help me write a new code. Thank you in advance.

Answer 1

[Sam Chak]

Hi @Elina

You likely want to find the value of x where the maximum value of occurs, but you may have inadvertently used the Newton–Raphson method to find the root where .

If you want to find the maximum point, you should apply the Newton–Raphson method to find the root of . By the way, the derivative of the function (f_prime) in your original code appears to be incorrect. Note that the symbolic result for the derivative is 100% correct! Please double check.

The f function should be and f_prime should be .

p           = @(x) 100*sqrt(x).*0.99.^(sqrt(x)) - 0.5*x;
dp          = @(x) 50*(0.99.^sqrt(x))./sqrt(x) + log(0.99)*50*0.99.^sqrt(x) - 0.5;
x           = linspace(1, 10000, 100000);

%% find the maximum point of p(x)
[val, idx]  = max(p(x))

val = 1.8759e+03

idx = 15994

xmax        = x(idx)

xmax = 1.6002e+03

%% find the root of p(x)
proot       = fzero(@(x) p(x), 7000)

proot = 7.2358e+03

%% Can use Newton–Raphson method to find the root of dp/dx
dproot      = fzero(@(x) dp(x), 1500)

dproot = 1.6001e+03

yyaxis left
plot(x, p(x))
grid on
yline(val,   '--')
xline(xmax,  '--')
xline(proot, '--')
xlabel({'$x$'}, 'interpreter', 'latex')
ylabel({'$p(x)$'}, 'interpreter', 'latex')

yyaxis right
plot(x, dp(x)), ylim([-0.75, 1]), grid on
ylabel({'$p''(x)$'}, 'interpreter', 'latex')

[Elina]

Hi @Sam Chak!

Thank you for the quick reply, and you are right, I was trying to find the root of derivative. And did not realise I started calculating the root of the function. But I don't realy understand the code you gave me, even thogh it seems to work.

I would really apriciate it if you could break down the code, that does everything do?

I have tride to add my understanding. I mostly don't understand

xroot = fzero(@(x) p(x), 7000

p           = @(x) 100*sqrt(x).*0.99.^(sqrt(x)) - 0.5*x; % Defines the funtion p(x)
dp          = @(x) 50*(0.99.^sqrt(x))./sqrt(x) - 0.502517*0.99.^sqrt(x) - 0.5; % Defines the funtion p'(x)
x           = linspace(1, 10000, 100000); % Kind of understan, creat a vector with 100,000 linearly spaced values ranging from 1 to 10,000.
dpdx        = gradient(p(x))./gradient(x); % Gives an approximation of the derivative p'(x) at each point. 
[val, idx]  = max(p(x)) %finds the maximum value of the function p(x)
xmax        = x(idx) % Where xmax is the max of the function p(x)
xroot       = fzero(@(x) p(x), 7000) % Result: the root of the derivative? How does it do that, it looks like it is finding the root of p(x).

Also how did you figure out the derivative? Did you calculate it yourself or did you compute it somhow?

But to find the root of the derivative the above code (in this comment) is all I need right? xroot is the root of the derivative?

I think what makes me confused is I don't need a loop and I don't recognize the formulas. Like what would gradient(p(x))./gradient(x) be in "paper and pen" math.

And how did you guess x0=7000, I dont really know where to start to gues the root of the function.

Anyway, appricate all the help I can get. I realize ther are a lot of qestions, I appricate it if you even just answer a few. I feel like if I understand every comand I will be able to remeber it better and re-use it. :)

[Sam Chak]

Hi @Elina

Since this is a numerical mathematics course and your Newton–Raphson code is capable of finding the root, there is nothing to fix. My explanation and demonstration using the built-in MATLAB functions (max and fzero) were intended to illustrate that you should find the root of ​ instead of .

The max() function is used to find both the maximum value of the graph and its indexed position. Note that indexing is also not covered in some numerical textbooks. The fzero() function is designed to find the root of a given function with an initial guess.

Regarding how I guessed the initial estimate of 7,000, which is relatively close to the true solution at , I simply plotted the graph of until I observed that the crossing at occurs between 6,000 and 8,000. Thus, I took the midpoint. This logical approach to determining the initial guess graphically is often taught by professors who emphasize logical reasoning, but it is not typically covered in most textbooks, nor is it mentioned on Wikipedia.

[Sam Chak]

As for how I determined the derivative of , you actually computed it correctly on your own. I simply took your result and created a function handle for dp, which stands for the derivative of .

dp = @(x) 50*(0.99.^sqrt(x))./sqrt(x) + log(0.99)*50*0.99.^sqrt(x) - 0.5;

Your code:

syms x
f   = inline('100 .* sqrt(x) .* 0.99.^(sqrt(x)) - 0.5 .* x', 'x')

f =

     Inline function:
     f(x) = 100 .* sqrt(x) .* 0.99.^(sqrt(x)) - 0.5 .* x

df  = diff(f(x), x)

[Walter Roberson]

It is recommended that you avoid using inline(). Replace it with

syms x
f   = 100 .* sqrt(x) .* 0.99.^(sqrt(x)) - 0.5 .* x;

df  = diff(f(x), x)