Making correct Jacobian for fsolve()

Hello,

I'm currently trying to use fSolve to solve a (relatively) simple nonlinear equation to represent the effects of drag on an airplane. The equation is the integral of x/(A*x^3 + B), where A and B are constants. I used wolfram to find the integral, and used it as the objective in fsolve. It worked just fine, but I wanted to try to speed it up, as I have to solve the equation repeatedly under a large number of conditions.

For this reason, I was going to use a jacobian to speed up calculations. I used the jacobian J = x/(A*x^3 + B). I was under the impression that since this equation has only one variable that would be all I needed to do. However, after implementing this jacobian fsolve stopped working entirely.

Setup for fsolve:

options = optimoptions('fsolve','Jacobian','on');
[finalSpeed fval] = fsolve(@(finalSpeed)straightLineDrag(distance,currentSpeed,finalSpeed,A,B),finalSpeedEstimate,options);

My function for F and J:

[F,J] = straightLineDrag(dist,V_0,V_f,A,B)
signA = 1;
signB = 1;
if A < 0
    signA = -1;
end
if B < 0
    signB = -1;
end

A_1_3 = signA*norm(A^(1/3));
A_2_3 = norm(A^(2/3));
B_1_3 = signB*norm(B^(1/3));
B_2_3 = norm(B^(2/3));

term_1_1 = -2*sqrt(3)*atan((1-(2*A_1_3*V_f)/B_1_3)/sqrt(3));
term_1_2 = -2*sqrt(3)*atan((1-(2*A_1_3*V_0)/B_1_3)/sqrt(3));

term_2_1 = -2*log(B_1_3 + A_1_3*V_f);
term_2_2 = -2*log(B_1_3 + A_1_3*V_0);

term_3_1 = log(B_2_3 - A_1_3*B_1_3*V_f + A_2_3*(V_f^2));
term_3_2 = log(B_2_3 - A_1_3*B_1_3*V_0 + A_2_3*(V_0^2));

term_4 = dist*(6*A_2_3*B_1_3);

F = term_1_1 - term_1_2 + term_2_1 - term_2_2 + term_3_1 - term_3_2 + term_4;

J = V_f/(A*(V_f^3) + B);
end

A and B values: (I don't think it matters, but I'm including it anyways)

A_prime = rho*A_S/(2*m);
A = (C_D_0 + C_D_alpha*abs(alpha))*A_prime;
B = g*sin(alpha);

Where everything is a constant except for alpha, which is determined from the angle of a straight line path (varies +/- 15 degrees)

[Walter Roberson]

That integral has all kinds of conditions, and the information you give isn't really enough for us to determine whether those conditions are being violated. In particular we can't tell whether A*x^3 + B might become 0 for some x.

Do you get an error message with fsolve() ?

[Torsten]

F = (term_1_1 - term_1_2 + term_2_1 - term_2_2 + term_3_1 - term_3_2)/term_4;
J = V_f/(A*(V_f^3) + B);

Best wishes

Torsten.

Answer 1

[Matt J]

Since you have only 1 unknown, it might help to use fzero instead, which requires no Jacobian calculations.