MATLAB Answers

Troubles with interaction problem

1 view (last 30 days)
raquel
raquel on 6 Dec 2013
Edited: raquel on 7 Dec 2013
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 Comments

bym
bym on 6 Dec 2013
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 on 7 Dec 2013
It is like the way you wrote last on my code, the way I pasted it here messed it up :)

Sign in to comment.

Accepted Answer

bym
bym on 6 Dec 2013
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 Comment

raquel
raquel on 7 Dec 2013
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?

Sign in to comment.

More Answers (0)

Sign in to answer this question.