vpasolve for three nonlinear equations inside a for loop
4 ビュー (過去 30 日間)
古いコメントを表示
I am trying to solve the system of the three non-linear equations by using vpasolve for different values of parameter U that varies from 0-0.5, while other parameters are fixed. My code is just showing output for U=0, then it's showing an error(Second argument must be a vector of symbolic variables).
clc
clear all;
syms x y w z eb
t = 0.2./pi;
d = 0.2./pi;
U = 0;
while (U<0.5)
e=U./2;
a = U./(pi.*t);
f = imag((U./d).*((((t./d)./(sqrt(1-z.^2)))+w)./((z.^2.*(1+2.*(t./d)./sqrt(1-z.^2)))-(2.*(t./d).*w./...
sqrt(1-z.^2))-(t./d).^2-(( -e + U.*x)./d).^2-w.^2)));
g = imag((z+(( -e + U.*x)./d)+((t./d).*z)./(sqrt(1-z.^2)))./((z.^2.*(1+2.*(t./d)./sqrt(1-z.^2)))-(2.*(t./d).*w./...
sqrt(1-z.^2))-(t./d).^2-(( -e + U.*x)./d).^2-w.^2));
h = imag((z+((-e + U.*y)./d)+((t./d).*z)./(sqrt(1-z.^2)))./((z.^2.*(1+2.*(t./d)./sqrt(1-z.^2)))-(2.*(t./d).*w./...
sqrt(1-z.^2))-(t./d).^2-((-e + U.*y)./d).^2-w.^2));
s=int(f,z,-Inf,0);
u=int(g,z,-Inf,0);
v=int(h,z,-Inf,0);
eq1=w-(-1./pi).*s==0;
eq2=y-(-1./pi).*u==0;
eq3=x-(-1./pi).*v==0;
% sol = vpasolve(eqs,vars);
[x,y,w] = vpasolve([eq1, eq2, eq3],[x,y,w],[1 0 0]);
%[sol.x sol.y sol.w]
%solutions = [solx,soly,solw]
n_up = double(x);
n_down = double(y);
d_ind = double(w);
m = abs(n_up-n_down);
eqn=((eb.*(1+2.*(t./d)./sqrt(1-eb)))-(2.*(t./d).*d_ind./sqrt(1-eb))-(t./d).^2-(( -e + U.*n_up)./d).^2-d_ind.^2);
e_abs = vpasolve(eqn,eb);
e_abs1=sqrt(e_abs);
e_abs2=-sqrt(e_abs);
weight = ((1./2).*(1-(e_abs1).^2).*(((sqrt(1-e_abs1.^2)).*(1+((( -e + U.*n_up)./d)./(e_abs1))))./(((1-(e_abs1).^2)...
.*((sqrt(1-e_abs1.^2))+2.*(t./d)))+((t./d).*e_abs1.^2)+(t./d).*(d_ind./d))));
fprintf('%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f\n', [a,n_up,n_down,m,d_ind,e_abs1,e_abs2,weight]');
if (U==0.5)
break
end
U = U+0.1;
end
採用された回答
Ameer Hamza
2020 年 6 月 14 日
In this line
[x,y,w] = vpasolve([eq1, eq2, eq3],[x,y,w],[1 0 0]);
you are overwriting the values of x, y, and w and converting them from symbolic to numeric. Use different variable names. See the following code
clc
clear all;
syms x y w z eb
t = 0.2./pi;
d = 0.2./pi;
U = 0;
while (U<0.5)
e=U./2;
a = U./(pi.*t);
f = imag((U./d).*((((t./d)./(sqrt(1-z.^2)))+w)./((z.^2.*(1+2.*(t./d)./sqrt(1-z.^2)))-(2.*(t./d).*w./...
sqrt(1-z.^2))-(t./d).^2-(( -e + U.*x)./d).^2-w.^2)));
g = imag((z+(( -e + U.*x)./d)+((t./d).*z)./(sqrt(1-z.^2)))./((z.^2.*(1+2.*(t./d)./sqrt(1-z.^2)))-(2.*(t./d).*w./...
sqrt(1-z.^2))-(t./d).^2-(( -e + U.*x)./d).^2-w.^2));
h = imag((z+((-e + U.*y)./d)+((t./d).*z)./(sqrt(1-z.^2)))./((z.^2.*(1+2.*(t./d)./sqrt(1-z.^2)))-(2.*(t./d).*w./...
sqrt(1-z.^2))-(t./d).^2-((-e + U.*y)./d).^2-w.^2));
s=int(f,z,-Inf,0);
u=int(g,z,-Inf,0);
v=int(h,z,-Inf,0);
eq1=w-(-1./pi).*s==0;
eq2=y-(-1./pi).*u==0;
eq3=x-(-1./pi).*v==0;
% sol = vpasolve(eqs,vars);
[xv,yv,wv] = vpasolve([eq1, eq2, eq3],[x,y,w],[1 0 0]);
%[sol.x sol.y sol.w]
%solutions = [solx,soly,solw]
n_up = double(xv);
n_down = double(yv);
d_ind = double(wv);
m = abs(n_up-n_down);
eqn=((eb.*(1+2.*(t./d)./sqrt(1-eb)))-(2.*(t./d).*d_ind./sqrt(1-eb))-(t./d).^2-(( -e + U.*n_up)./d).^2-d_ind.^2);
e_abs = vpasolve(eqn,eb);
e_abs1=sqrt(e_abs);
e_abs2=-sqrt(e_abs);
weight = ((1./2).*(1-(e_abs1).^2).*(((sqrt(1-e_abs1.^2)).*(1+((( -e + U.*n_up)./d)./(e_abs1))))./(((1-(e_abs1).^2)...
.*((sqrt(1-e_abs1.^2))+2.*(t./d)))+((t./d).*e_abs1.^2)+(t./d).*(d_ind./d))));
fprintf('%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f\n', [a,n_up,n_down,m,d_ind,e_abs1,e_abs2,weight]');
if (U==0.5)
break
end
U = U+0.1;
end
2 件のコメント
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Symbolic Math Toolbox についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!