Troubles with interaction problem

조회 수: 5 (최근 30일)
raquel
raquel 2013년 12월 6일
편집: raquel 2013년 12월 7일
Hey guys! I'm having trouble finding a value for "beta" in the following code:
e=1.239e-7;
c=5.154e-4;
d=4.823;
me=2.0475e8;
mv=28.96;
se=0.06.*me;
sv=6.081;
imax=15; %Max number of interactions
tol=0.001; %Tolerance for convergence
Edes(1)=me;
Vdes(1)=mv;
Ered_des(1)=(Edes(1)-me)./se;
Vred_des(1)=(Vdes(1)-me)./se;
dg_de(1)=se./(me+se.*Ered_des(1));
dg_dv(1)=-2.*c.*sv.*(mv+sv.*Vred_des(1))./(d+c.*(mv+sv.*Vred_des(1)).^2);
alfa_e(1)=(dg_de(1))./(sqrt(dg_de(1).^2+dg_dv(1).^2));
alfa_v(1)=(dg_dv(1))./(sqrt(dg_de(1).^2+dg_dv(1).^2));
syms beta
g(1)=log( e.*((-alfa_e(1).*beta(1)).*se+me)./(c.*(-alfa_v(1).*beta(1).*sv+mv).^2+d) );
[beta(1)] = solve(g(1)==0);
beta(1)=vpa((subs(beta(1))));
for ii=1:imax
beta(ii+1)=0;
if abs(beta(ii+1)-beta(ii))>0.001
Ered_des(ii+1)=-alfa_e(ii).*beta(ii);
Vred_des(ii+1)=-alfa_v(ii).*beta(ii);
dg_de(ii+1)=se./(me+se.*Ered_des(ii+1));
dg_dv(ii+1)=-2.*c.*sv.*(mv+sv.*Vred_des(ii+1))./(d+c.*(mv+sv.*Vred_des(ii+1)).^2);
alfa_e(ii+1)=(dg_de(ii+1))./(sqrt(dg_de(ii+1).^2+dg_dv(ii+1).^2));
alfa_v(ii+1)=(dg_dv(ii+1))./(sqrt(dg_de(ii+1).^2+dg_dv(ii+1).^2));
syms beta
g(ii+1)=log( e.*((-alfa_e(ii+1).*beta(ii+1)).*se+me)./(c.*(-alfa_v(ii+1).*beta(ii+1).*sv+mv).^2+d) );
[beta(ii+1)] = solve(g(ii+1)==0);
beta(ii+1)=vpa((subs(beta(ii+1))));
else
Beta=double(beta(ii+1));
end
if ii==imax
disp('Exceeded max. number of interactions!')
end
end
When I try this, the following error appears:
"Error using mupadmex Error in MuPAD command: Index exceeds matrix dimensions.
Error in sym/subsref (line 1577) B = mupadmex('symobj::subsref',A.s,inds{:});
Error in InvarianteHL (line 62) g(ii+1)=log( e.*((-alfa_e(ii+1).*beta(ii+1)).*se+me)./(c.*(-alfa_v(ii+1).*beta(ii+1).*sv+mv).^2+d) );"
Please please help?
  댓글 수: 2
bym
bym 2013년 12월 6일
편집: bym 2013년 12월 6일
this line
for ii=1:imax beta(ii+1)=0;
is improper syntax maybe there is a missing return and you intend
for ii=1:imax
beta(ii+1)=0;
raquel
raquel 2013년 12월 7일
편집: raquel 2013년 12월 7일
It is like the way you wrote last on my code, the way I pasted it here messed it up :)

댓글을 달려면 로그인하십시오.

채택된 답변

bym
bym 2013년 12월 6일
here is a partial answer, without using symbolic tool box...see my comment above and you should be able to adapt the code to what you want
clear
e=1.239e-7; c=5.154e-4; d=4.823;
me=2.0475e8; mv=28.96;
se=0.06.*me; sv=6.081;
imax=15; %Max number of interactions
tol=0.001; %Tolerance for convergence
Edes(1)=me;
Vdes(1)=mv;
Ered_des(1)=(Edes(1)-me)./se;
Vred_des(1)=(Vdes(1)-me)./se;
dg_de(1)=se./(me+se.*Ered_des(1));
dg_dv(1)=-2.*c.*sv.*(mv+sv.*Vred_des(1))./(d+c.*(mv+sv.*Vred_des(1)).^2);
alfa_e(1)=(dg_de(1))./(sqrt(dg_de(1).^2+dg_dv(1).^2));
alfa_v(1)=(dg_dv(1))./(sqrt(dg_de(1).^2+dg_dv(1).^2));
%syms beta
g=@(beta)log( e.*((-alfa_e(1).*beta).*se+me)./(c.*(-alfa_v(1).*beta.*sv+mv).^2+d) );
%[beta(1)] = solve(g(1)==0);
%beta(1)=vpa((subs(beta(1))));
b(1) = fzero(g,10);
  댓글 수: 1
raquel
raquel 2013년 12월 7일
편집: raquel 2013년 12월 7일
From the line:
syms beta
g(ii+1)..
I replaced with:
g=@(beta)log( e.*((-alfa_e(ii+1).*beta).*se+me)./(c.*(-alfa_v(ii+1).*beta.*sv+mv).^2+d) );
beta(ii+1) = fzero(g,10);
dif=abs(beta(ii+1)-beta(ii));
if dif<0.001
Beta=double(beta(ii+1))
break
end
And it works! :D Thank you.
But in this case I'd have to make a guess for beta, right? As you did using fzero(g,10).. Isn't there any other way to solve for beta?

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Conversion Between Symbolic and Numeric에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by