vpasolve for three nonlinear equations inside a for loop
2 views (last 30 days)
Show older comments
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
Accepted Answer
Ameer Hamza
on 14 Jun 2020
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 Comments
More Answers (0)
See Also
Categories
Find more on Symbolic Math Toolbox 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!