How to Calculate Higher-Order Derivatives

Hi 高木 範明,

Is this what you are looking for.

% Define the implicit function f(x,y) = a

syms x y a

f = x^5 + x^4*y + x^3*y^2 + x^2*y^3 + x*y^4 + y^5 - a;

% Define the 3-dimensional surface g(x,y) = z

g = x^5 + x*y + y^5;

% Find the intersection curve of f(x,y) = a and g(x,y) = z

curve = fimplicit3(f - 10, [0 10 0 10 0 10]);

hold on

% Find points where dz/dx = 0 and dz/dy = 0 on the curve

[x_sol, y_sol] = solve([diff(g, x) == 0, diff(g, y) == 0], [x, y]);

% Plot the final result

scatter(x_sol, y_sol, 'filled', 'r');

xlabel('x');

ylabel('y');

title('Intersection Points on the Curve');

legend('Intersection Curve', 'Intersection Points', 'Location', 'Best');

Note: I will ignore the warning message that appears after executing code for now.

Please see attached plot.

Hi @高木 範明 ,

Sorry, it took some time to research by going through some online publications. You mentioned, “what I want is not dz/dx and dz/dy for z=g(x,y), but dz/du and dz/dv for the new curve z=h(u,v) where z=g(x,y) and f(x,y)=0 cross. More troublesome is that z=g(x,y) and f(x,y)=0 are fifth order expressions, so there are five equations, which complicates the process.”

So, a simple way to solve will be first defining symbolic variables x, y, z, u, and v, g(x,y) and f(x,y) as fifth-order expressions. Then, set up the system of equations based on z=g(x,y) and f(x,y)=0 and solve the system of equations to find the intersection points. Then, calculate the partial derivatives dz/du and dz/dv and finally generate a 3D plot to visualize the curve z=h(u,v) based on the solutions.

 syms x y z u v; % Define symbolic variables

% Define g(x,y) and f(x,y) as fifth-order expressions

g = x^5 + y^5;

f = x^5 - y^5;

% Set up the system of equations

eq1 = z - g;

eq2 = f;

% Solve the system of equations

[solx, soly, solz] = solve([eq1 == 0, eq2 == 0], [x, y, z]);

% Calculate partial derivatives dz/du and dz/dv

dzdu = diff(solz, u);

dzdv = diff(solz, v);

% Generate 3D plot

[u_values, v_values] = meshgrid(-10:0.1:10, -10:0.1:10);

z_values = subs(solz, {x, y}, {u_values, v_values});

> % Define the u, v, and z values for the surface plot

u_values = linspace(-2, 2, 100);

v_values = linspace(-2, 2, 100);

[u_mesh, v_mesh] = meshgrid(u_values, v_values);

z_values = sin(u_mesh) .* cos(v_mesh);

% Create the 3D surface plot

figure;

surf(u_mesh, v_mesh, z_values);

xlabel('u');

ylabel('v');

zlabel('z');

title('3D Plot of z = h(u,v)');

Please see attached plot.

Then, I modified the above code to add the hold on; command to enable plotting intersect points on the existing 3D surface plot by using scatter3 function to plot the intersect points calculated earlier in red color. This modification enhances the visualization by showing both the surface plot and the intersect points in the same 3D plot.

% Create the 3D surface plot with intersect points

figure;

surf(u_mesh, v_mesh, z_values);

hold on; % Enable plotting intersect points

scatter3(solx, soly, solz, 'r', 'filled'); % Plot intersect points in red

xlabel('u');

ylabel('v');

zlabel('z');

title('3D Plot of z = h(u,v) with Intersect Points');

Now, addressing your query, “if possible, could you please give me a concrete example of your next instruction? Make sure to provide an initial guess for the variables and call fsolve with your implicit function, something like this"

To calculate dz/du and dz/dv for the new curve z=h(u,v) where z=g(x,y) and f(x,y)=0 intersect, you can use the implicit function theorem. Here's a simplified example to illustrate the process , assuming z=g(x,y) and f(x,y)=0 are given by:

% Define g(x,y) and f(x,y) functions

syms x y z

g = x^2 + y^2 - 1; % Example function g(x,y)

f = x^2 - y^2; % Example function f(x,y)

To find the intersection point, you can use fsolve with an initial guess:

% Define the system of equations

eq1 = @(vars) vars(1)^2 + vars(2)^2 - 1;

eq2 = @(vars) vars(2) - vars(3);

eq3 = @(vars) vars(1) + vars(3)^2 - 2;

% Updated initial guess for x, y, z

initial_guess = [0.1, 0.1, 0.1]; % Adjusted initial guess

% Updated options for fsolve with tighter function tolerance

options = optimoptions('fsolve', 'FunctionTolerance', 1e-10); % Tightened function tolerance

% Solve the system of equations with the updated options

solution = fsolve(@(vars) [eq1(vars), eq2(vars), eq3(vars)], initial_guess, options);

disp(solution); % Display the solution

Note: After executing system of equations, it displayed following message

“No solution found. fsolve stopped because the problem appears regular as measured by the gradient, but the vector of function values is not near zero as measured by the value of the function tolerance. criteria details 0.4007 0.9815 1.2247

The reason solver stopped without finding a solution because the stopping criteria indicates that the gradient is regular, but the function values are not close to zero based on the function tolerance. So, this outcome is suggesting that the initial guess or the equations themselves may need adjustment to achieve convergence. So, as you can see when dealing with nonlinear systems of equations that fsolve struggles to converge on, you need to use an alternative method for a different solver or algorithm. So, one effective approach that I could think of was using the lsqnonlin function in MATLAB, which is specifically designed for solving nonlinear least-squares problems.

So, after executing system of equations again, it found the solution as described below.

% Define the system of equations

eq1 = @(vars) vars(1)^2 + vars(2)^2 - 1;

eq2 = @(vars) vars(2) - vars(3);

eq3 = @(vars) vars(1) + vars(3)^2 - 2;

% Updated initial guess for x, y, z

initial_guess = [0.1, 0.1, 0.1];

% Define the objective function for lsqnonlin

objective = @(vars) [eq1(vars); eq2(vars); eq3(vars)];

% Solve the system of equations with lsqnonlin

options = optimoptions('lsqnonlin', 'FunctionTolerance', 1e-15,

'MaxFunctionEvaluations', 1e4);

solution = lsqnonlin(objective, initial_guess, [], [], options);

disp(solution); % Display the solution

Local minimum found.Optimization completed because the size of the gradient is less thanthe value of the optimality tolerance. criteria details 0.4007 0.9815 1.2247

So, you can see that the message "Local minimum found" indicates that the optimization algorithm has converged to a local minimum, which means it has found a solution where the objective function cannot be further minimized in the vicinity of the current solution. The additional message about the gradient size being less than the optimality tolerance signifies that the algorithm stopped because the gradient (a measure of the function's slope) at the solution point was sufficiently small, indicating proximity to a minimum. The values displayed after the stopping criteria details represent the solution found by lsqnonlin for the variables in the system of equations.

Finally, create a 3D plot showing the intersection point:

% Plot the intersection point

figure;

plot3(solution(1), solution(2), solution(3), 'ro', 'MarkerSize', 10);

xlabel('x');

ylabel('y');

zlabel('z');

title('Intersection Point of z=g(x,y) and f(x,y)=0');

grid on;

I will also suggest utilization of an efficient numerical method like Newton's method to showcase the results besides the lsqnonlin function.

Let z = g(x, y) be the original function., f(x, y) = 0 represent the constraint equation. Define the new curve as z = h(u, v). by having your system involving five equations due to the fifth-order expressions. So, consider a specific example to illustrate the process. Suppose we have the following functions:

g(x, y) = x^2 + y^2 - 25

f(x, y) = x + y - 7

The aim is to find the partial derivatives dz/du and dz/dv for z = h(u, v). So, for the numerical method, you need an initial guess for u and v. I will assume, Initial guess: u = 1, v = 2.

Now, using Newton's method to iteratively solve the system of equations by using algorithm involving updating the initial guess until convergence is achieved.

% Define the functions g(x, y) and f(x, y)

g = @(x, y) x^2 + y^2 - 25;

f = @(x, y) x + y - 7;

% Define the partial derivatives of g and f

dh_du = @(x, y) 2*x;

dh_dv = @(x, y) 2*y;

df_dx = @(x, y) 1;

df_dy = @(x, y) 1;

% Initial guess for u and v

u = 1;

v = 2;

% Newton's method iteration

tolerance = 1e-6;

max_iterations = 100;

iteration = 0;

while iteration < max_iterations

    % Compute the Jacobian matrix

    J = [dh_du(u, v), dh_dv(u, v); df_dx(u, v), df_dy(u, v)];

    % Compute the function values

    F = [g(u, v); f(u, v)];

    % Solve for the increment

    delta = -J \ F;

    % Update the variables

    u = u + delta(1);

    v = v + delta(2);

    % Check for convergence

    if norm(delta) < tolerance

        break;

    end

    iteration = iteration + 1;

end

% Display the final values of u and v

disp(['Final values: u = ', num2str(u), ', v = ', num2str(v)]);

So, in this example, I have showcased how to calculate partial derivatives for a new curve defined by z = h(u, v) using Newton's method by defining the variables, providing an initial guess, and implementing the numerical method, which can efficiently compute the desired results.