Solving an implicit equation with fsolve and validating with fplot

10 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.

Categories

Find more on Mathematics in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!