Solving an implicit equation with fsolve and validating with fplot

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)])

 採用された回答

Ryan McBurney
Ryan McBurney 2021 年 2 月 1 日

0 投票

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.

その他の回答 (1 件)

Alan Weiss
Alan Weiss 2021 年 1 月 31 日
編集済み: Alan Weiss 2021 年 1 月 31 日

0 投票

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 件のコメント

Ryan McBurney
Ryan McBurney 2021 年 2 月 1 日
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])

サインインしてコメントする。

カテゴリ

ヘルプ センター および File ExchangeCreating and Concatenating Matrices についてさらに検索

製品

Community Treasure Hunt

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

Start Hunting!

Translated by