Solving an implicit equation with fsolve and validating with fplot

20 views (last 30 days)
Below is my code to try and solve this function:
D/(log10(RX*cfe)-C))^2 = cfe
I'm not sure if I'm using fsolve correctly. I get a value for cfe as 0.0016 but when I try to some validation by plotting two sides of the equation the curves cross at .0017377. That's about a 7% difference in values so I am questioning if I'm using fsove correctly.
I'm actually unsure of how to implement the initial guess which I call cfi.
I thought maybe adjusting the options would make the code run longer and would get a closer answer to the plotting method, but I'm not even sure if I'm doing that right.
Any direction would be helpful.
% AL, BE, E, and TWTD are all constants, which makes D and C also constants.
D = 0.242*(asin(AL)+asin(BE))/sqrt(E);
C = 1.26*log10(TWTD);
RX = 3.5e6 % Constant
options = optimoptions('fsolve','FunctionTolerance',1.0000e-1);
cfi = .455/(log10(RX))^2.58; % Initial guess
f = @(cfe,cfi) (D/(log10(RX*cfi)-C))^2 - cfe;
fsol = fsolve(@(cfe) f(cfe,cfi),0,options);
%------------------------------------------------------
syms cfe
eqnLeft = 0.242*(asin(AL)+asin(BE))/(log10(RX*cfe)-C);
eqnRight = sqrt(cfe)*sqrt(E);
fplot([eqnLeft eqnRight])
title([texlabel(eqnLeft) ' = ' texlabel(eqnRight)])

Accepted Answer

Ryan McBurney
Ryan McBurney on 1 Feb 2021
Fixed my problem with the following adjustments to my code.
cfi = .455/(log10(RX))^2.58;
f = @(cfe,RX,Me) ((D)/(log10(RX*cfe)-C))^2 - cfe;
fsol = fsolve(@(cfe) f(cfe,RX,Me),cfi);
I had a suspicion I was using the initial guess the wrong way and I believe that is what is giving me a correct answer. The answer now matches the fplot result.

More Answers (1)

Alan Weiss
Alan Weiss on 31 Jan 2021
Edited: Alan Weiss on 31 Jan 2021
You have set an enormous value of the FunctionTolerance option. Don't set any options and I expect that you will get a better answer.
If this advice does not help, then please give numeric values for your constants so that we might have a chance at running your code and seeing what is happening.
Alan Weiss
MATLAB mathematical toolbox documentation
  1 Comment
Ryan McBurney
Ryan McBurney on 1 Feb 2021
Here is the full code. Eventually I want to make Me and REN a vector of inputs.
clear all; clc;
Me = 3.5;
REN = 1.3e6;
LCFT = 12;
g = 1.4;
RX = REN*LCFT;
E = (g-1)/2 * Me^2;
TWTD = 1 + 0.89 * E;
A = E/TWTD;
B = (1+E)/(TWTD) - 1;
AL = (2*A-B)/sqrt(B^2+4*A);
BE = B/sqrt(B^2+4*A);
D = 0.242*(asin(AL)+asin(BE))/sqrt(E);
C = 1.26*log10(TWTD);
cfi = .455/(log10(RX))^2.58; % Initial guess
f = @(cfe,cfi) (D/(log10(RX*cfi)-C))^2 - cfe;
fsol = fsolve(@(cfe) f(cfe,cfi),0);
%-------------------------------------------
syms cfe
eqnLeft = 1/(log10(RX*cfe)-C);
eqnRight = sqrt(cfe)*sqrt(E)/(0.242*(asin(AL)+asin(BE)));
fplot([eqnLeft eqnRight])

Sign in to comment.

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!