Is there any problem with secant method in this code as i am not getting the required plot

31 Jan 2023
1 Réponse
Mise à jour 31 Jan 2023
1 Vue (30 jours)
Ran in:

0 votes

Ran in:
clear
i=1;
error=0.0001;
C=zeros(100,1);
aC=zeros(100,1);
for a=0.01:0.01:1
span=3.6;
yspan=[-span span];
k=0;
u0=[0.09,a*0.09];
c0=complex(0,0.1);
sol=ode45(@(y,u) func(y,u,c0,a),yspan,u0);
U=deval(sol,span);
A0=(U(2)/U(1))+a;
c1=complex(0,0.12);
sol=ode45(@(y,u) func(y,u,c1,a),yspan,u0);
U=deval(sol,span);
A1=(U(2)/U(1))+a;
c=(((c0*A1)-(c1*A0)/(A1-A0)));
while k<1000
sol=ode45(@(y,u) func(y,u,c,a),yspan,u0);
U=deval(sol,span);
A0=A1;
A1=(U(2)/U(1))+a;
c0=c1;
c1=c;
errR=abs(real(c1-c0));
errI=abs(imag(c1-c0));
if((errR<error)&&(errI<error))
break;
end
c=(((c0*A1)-(c1*A0)/(A1-A0)));
k=k+1;
end
c;
k
C(i)=imag(c);
aC(i)=a*imag(c);
i=i+1;
end
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
k = 1000
t=linspace(0.01,1,100)';
plot(t,C)
hold on
plot(t,aC)
hold off
ylim([0 1.1]);
xlim([0.01 1]);
legend('Ci','αCi');
function dudy=func(y,u,c,a)
dudy=zeros(2,1);
dudy(1,1)=u(2);
dudy(2,1)=(a*a+2*tanh(y)*(tanh(y)*tanh(y)-1)/(tanh(y)-c))*u(1);
end

Connectez-vous pour répondre à cette question.

Réponses (1)

Walter Roberson
Walter Roberson le 31 Jan 2023
Ran in:

0 votes

Ran in:
Your odes are generating infinite results early on. You are getting NaN for all c values.
WIth the below small modification to break early when NaN is detected, infinite values show up for c
isfirstnan = true;
i=1;
error=0.0001;
C=zeros(100,1);
aC=zeros(100,1);
for a=0.01:0.01:1
span=3.6;
yspan=[-span span];
k=0;
u0=[0.09,a*0.09];
c0=complex(0,0.1);
sol=ode45(@(y,u) func(y,u,c0,a),yspan,u0);
U=deval(sol,span);
A0=(U(2)/U(1))+a;
c1=complex(0,0.12);
sol=ode45(@(y,u) func(y,u,c1,a),yspan,u0);
U=deval(sol,span);
A1=(U(2)/U(1))+a;
c=(((c0*A1)-(c1*A0)/(A1-A0)));
while k<1000
sol=ode45(@(y,u) func(y,u,c,a),yspan,u0);
U=deval(sol,span);
A0=A1;
A1=(U(2)/U(1))+a;
c0=c1;
c1=c;
errR=abs(real(c1-c0));
errI=abs(imag(c1-c0));
if((errR<error)&&(errI<error))
break;
end
c=(((c0*A1)-(c1*A0)/(A1-A0)));
if (any(~isfinite(real(c))) || any(~isfinite(imag(c))))
if isfirstnan
fprintf('got first non-finite for c at a = %g, k = %g, i = %g\n', a, k, i);
disp(c);
isfirstnan = false;
end
break
end
k=k+1;
end
c;
k;
C(i)=imag(c);
aC(i)=a*imag(c);
i=i+1;
end
got first non-finite for c at a = 0.01, k = 13, i = 1
-Inf + Infi
whos C, min(C), max(C)
Name Size Bytes Class Attributes C 100x1 800 double
ans = -Inf
ans = Inf
whos aC, min(aC), max(aC)
Name Size Bytes Class Attributes aC 100x1 800 double
ans = -Inf
ans = Inf
%{
t=linspace(0.01,1,100)';
plot(t,C)
hold on
plot(t,aC)
hold off
ylim([0 1.1]);
xlim([0.01 1]);
legend('Ci','αCi');
%}
function dudy=func(y,u,c,a)
dudy=zeros(2,1);
dudy(1,1)=u(2);
dudy(2,1)=(a*a+2*tanh(y)*(tanh(y)*tanh(y)-1)/(tanh(y)-c))*u(1);
end

Catégories

En savoir plus sur PHY Components dans Centre d'aide et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by