vpasolve didn't return real value

2 次查看（过去 30 天）

Hadi Norouzi 2019-5-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/463390-vpasolve-didn-t-return-real-value

评论： Jan 2019-5-23

input.xlsx

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)