vpasolve didn't return real value
2 views (last 30 days)
Show older comments
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 Comments
Jan
on 23 May 2019
@Hadi Norouzi: Please comments in the section for comments. Accepting an answer means, that the problem is solved.
Answers (0)
See Also
Categories
Find more on Get Started with 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!