Filter löschen
Filter löschen

fzero coupled with ode45

1 Ansicht (letzte 30 Tage)
pritha
pritha am 23 Jan. 2023
Kommentiert: Walter Roberson am 27 Jan. 2023
I'm getting, the following comments while trying to solve ....
rspan = [0 10];
y0 =[0 0 0 0 0 0 0 0];
[r,y] = ode45(@(r,y) fun(r,y), rspan, y0)
Error using fzero
Function values at the interval endpoints must differ in sign.

Error in solution>fun (line 24)
Opp = fzero(@(Op)fp(Op), [-10 10]);

Error in solution (line 3)
[r,y] = ode45(@(r,y) fun(r,y), rspan, y0)

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
function F = fun(r,y)
%% Parameters
X1 = 100;
X2 = 100;
M = 151;
S = 92;
R = 0.8/197.3;
m1 = 110;
g1 = 6.26;
m2 = 120;
g2 = 8.17;
m3 = 130;
g3 = 9.33;
Q = .5;
a = 7.33*10^(-4);
X11 = X1 - g1*(1-(a*g1*y(1)/2))*y(1);
X21 = X2 - g1*(1-(a*g1*y(1)/2))*y(1);
%% Evaluation of fp and fn
fp = @(Op)(1 - (X11.*R) - sqrt(Op.^2 + (X11.*R).^2) - Op./tan(Op));
Opp = fzero(@(Op)fp(Op), [-10 10]);
fn = @(On)(1 - (X21.*R) - sqrt(On.^2 + (X21.*R).^2) - On./tan(On));
Onn = fzero(@(On)fn(On),[-10 10]);
Ep1 = sqrt((Opp/R).^2 + X11^2);
En1 = sqrt((Onn/R).^2 + X21^2);
%% Functions
Wp = sqrt(Ep1.^2 - X11^2);
Wn = sqrt(En1.^2 - X21^2);
ap = Wp.*r;
an = Wn.*r;
Jp = besselj(0,ap);
Jp1 = besselj(1,ap);
Jn = besselj(0,an);
Jn1 = besselj(1,an);
%% Component
Kp = sqrt((Ep1.^2 + X11^2)./Ep1.^2) .* Jp .* r;
Lp = sqrt((Ep1.^2 - X11^2)./Ep1.^2) .* Jp1 .* r;
Kn = sqrt((En1.^2 + X21^2)./En1.^2) .* Jn .* r;
Ln = sqrt((En1.^2 - X21^2)./En1.^2) .* Jn1 .* r;
%% Constant
Ap = Wp*R;
An = Wn*R;
NJp = besselj(0,Ap);
NJn = besselj(0,An);
Np = (S/(4*pi*R^3))^(1/2) .* (Ep1.*(Ep1-X11)*R).^(1/2)./(NJp.*(2*R.*Ep1.^2 - 2.*Ep1 + X11).^(1/2));
Nn = (M/(4*pi*R^3))^(1/2) .* (En1.*(En1-X21)*R).^(1/2)./(NJn.*(2*R.*En1.^2 - 2.*En1 + X21).^(1/2));
p1 = ((Np.^2 .* (Kp.^2 - Lp.^2)) + (Nn.^2 .* (Kn.^2 - Ln.^2)))./(4*pi*r.^2);
p2 = ((Np.^2 .* (Kp.^2 + Lp.^2)) + (Nn.^2 .* (Kn.^2 + Ln.^2)))./(4*pi*r.^2);
p3 = ((Np.^2 .* (Kp.^2 + Lp.^2)) - (Nn.^2 .* (Kn.^2 + Ln.^2)))./(4*pi*r.^2);
p4 = (Np.^2 .* (Kp.^2 + Lp.^2))./(4*pi*r.^2);
%% Equations
F(1) = y(2);
F(2) = - (2*y(2)/r) + (m1^2*y(1)) - (g1*(1- a*g1*y(1))*p1);
F(3) = y(4);
F(4) = - (2*y(4)/r) + (m2^2*y(3)) - (g2*p2);
F(5) = y(6);
F(6) = - (2*y(6)/r) + (m3^2*y(5)) - (g3*p3/2);
F(7) = y(8);
F(8) = - (2*y(8)/r) - (Q*p4);
F=[F(1);F(2);F(3);F(4);F(5);F(6);F(7);F(8)];
end
Error using fzero (line 274)
Function values at the interval endpoints must differ in sign.
Error in error_code>fun (line 26)
Opp = fzero(@(Op)fp(Op), [-10 10]);
Error in error_code>@(r,y)fun(r,y) (line 5)
[r,y] = ode45(@(r,y) fun(r,y), rspan, y0)
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in error_code (line 5)
[r,y] = ode45(@(r,y) fun(r,y), rspan, y0)
Warning: Function fields has the same name as a MATLAB builtin. We suggest you rename the function
to avoid a potential name conflict.
Please help me out to solve this problem, where fzero is coupled with ode45.
  8 Kommentare
pritha
pritha am 27 Jan. 2023
I'm sorry. It is X11. That means,
fp = vpasolve(((Op./((1 - (X11.*R) - sqrt(Op.^2 + (X11.*R).^2)).*tan(Op))))==tan(Op), Op).
Walter Roberson
Walter Roberson am 27 Jan. 2023
You still have a division by tan(Op) . That leads to an indefinite number of solutions that cannot be chosen between without additional information.
If you move out that tan(Op) to form
((Op./((1 - (X11.*R) - sqrt(Op.^2 + (X11.*R).^2))))) == tan(Op)^2
then by examination we can see that Op == 0 is a solution except in the case that X11 == 0 or R == 0.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Walter Roberson
Walter Roberson am 23 Jan. 2023
fp = @(Op)(1 - (X11.*R) - sqrt(Op.^2 + (X11.*R).^2) - Op./tan(Op));
Opp = fzero(@(Op)fp(Op), [-10 10]);
I was sure that I already did a write-up about this about two days ago, but I cannot find the relevant post.
You have a division by tan(Op) . That is going to be a problem if tan(Op) ever passes through zero. Which it does in the range [-10 10] at Op = -3*pi, -2*pi, -pi, 0, pi, 2*pi, and 3*pi .
fzero cannot handle functions that have a discontinuity in the search range.
When I looked at a previous post of yours (that I cannot find now), I suspected that you wanted to be working in degrees rather than radians, such as -10 degrees to +10 degrees. That would require using tand() instead of tan(). However it would still have a discontinuity at 0.
Yes it is true that limit(Op ./ tan(Op), op, 0) is 1 -- but MATLAB does not use limits in calculations unless you specifically ask for limits to be used symbolically.

Weitere Antworten (0)

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Tags

Produkte


Version

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by