vpasolve didn't return real value
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Hello. I have code that use vpasolve to solve equation, but after a short time that loops goes on the answer that vpa solve return is unreal( and i know it should be real).
is there any way to be sure only real value return in Variable? thanks
clc;clear all;close all;
%%initial value
B=xlsread('input','input','A2');
s0=xlsread('input','input','B2');counter=1;
delt=xlsread('input','input','C2');
delx=xlsread('input','input','D2');
teta=xlsread('input','input','E2');
a=xlsread('input','input','H2');
b=xlsread('input','input','I2');
T=xlsread('input','input','J2');
g=xlsread('input','input','K2');
X0=xlsread('input','input','B4:Q4');
Q0=xlsread('input','input','B5:Q5');
y=xlsread('input','input','B6:Q6');
%% main program
maty=zeros(200,16);
matQ=zeros(200,16);
matv=zeros(200,16);
for t=0.005:delt:T
for i = 1:16
v(i)=Q0(i)/(B*y(i));
end
for i = 1:15
syms x
[ys(i)]=vpasolve(0==(x-y(i))+(teta*(y(i+1)-y(i))/(1+teta*(v(i+1)-v(i)))*(v(i)-(g*x)^0.5)));
vs(i)=(v(i)-teta*((g*y(i+1))^0.5*v(i)-(g*y(i))^0.5*v(i+1)))/(1-teta*(v(i)-v(i+1)-(g*y(i))^0.5+(g*y(i+1))^0.5));
Qs(i)=B*ys(i)*vs(i);
Sfs(i)= a* vs(i)+ b * vs(i)^2;
[yr(i)]=vpasolve(0==(y(i+1)-x)-teta*(y(i+1)-y(i))/(1+teta*(v(i+1)-v(i)))*(v(i+1)+(g*x)^0.5));
vr(i)=(v(i+1)-teta*((g*y(i))^0.5*v(i+1)-(g*y(i+1))^0.5*v(i)))/(1+teta*(v(i+1)-v(i)+(g*y(i+1))^0.5-(g*y(i))^0.5));
Qr(i)=B*yr(i)*vr(i);
Sfr(i)=a* vr(i)+ b * vr(i)^2;
end
%upstream condition
Qp(1)= 0.0003+0.0075*exp(-0.5*abs((t-50.7261)/24.495)^2.9063);
syms X
[ypp]=vpasolve(0==((vs(1)-2*(g*ys(1))^0.5)+g*(s0-Sfs(1))*delt)+(2*(g*X)^0.5)-(Qp(1)/(B*X)));
for i=1:14
yppp(i) = 0.5*(yr(i)+ys(i+1))+ teta*((Qr(i)-Qs(i+1))/(2*B));
end
%downstream condition
ypppp=0.1201*exp(-0.5*((t-81.7891)/23.2503)^2);
%start velocity
vpp=Qp(1)/(B*ypp);
%end velocity
vpppp= (vr(15)+2*(g*yr(15))^0.5)+(g*(s0-Sfr(15))*delt)-2*(g*ypppp)^0.5;
for i = 1:14
vppp(i)=0.5*(vr(i)+vs(i+1))+0.5*teta*((vr(i)^2-vs(i+1)^2)+g*(yr(i)-ys(i+1)))+g*delt*(s0-0.5*(Sfr(i)+Sfs(i+1)));
end
%combine yp & vp
yp=horzcat(ypp,yppp,ypppp);
vp=horzcat(vpp,vppp,vpppp);
if t==0:1:200
plot(X0,yp)
end
for i=1:15
Qp(i+1)=B*vp(i+1)*yp(i+1);
end
%output
if t==(1*counter)
filename = 'y.xlsx';
filename1= 'V.xlsx';
filename2= 'Q.xlsx';
A=double(yp);
C=double(vp);
D=double(Qp);
for i=1:16
maty(counter,i)=A(1,i);
matQ(counter,i)=D(1,i);
matv(counter,i)=C(1,i);
my{counter,i}=[A(1,i)];
mv{counter,i}=[C(1,i)];
mq{counter,i}=[D(1,i)];
end
end
Q0=Qp;
y=yp;
aaaa=aaaa+1
end
xlswrite(filename,my)
xlswrite(filename1,mv)
xlswrite(filename2,mq)
2 comentarios
Jan
el 23 de Mayo de 2019
@Hadi Norouzi: Please comments in the section for comments. Accepting an answer means, that the problem is solved.
Respuestas (0)
Ver también
Categorías
Más información sobre Get Started with Symbolic Math Toolbox en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!