# Troubles with interaction problem

3 ビュー (過去 30 日間)
raquel 2013 年 12 月 6 日

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) );"
##### 2 件のコメント表示 1 件の古いコメント非表示 1 件の古いコメント
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 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 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?

サインインしてコメントする。

### カテゴリ

Find more on Common Operations in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by