Troubles with interaction problem

1 次查看（过去 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 Center 和 File Exchange 中查找有关 Common Operations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by