DOUBLE cannot convert the input expression into a double array.

Hi, I wrote the following code to generate the smmoth curves. When executed, the program shows the following errors:

The following error occurred converting from sym to double:

Error using symengine (line 58)

DOUBLE cannot convert the input expression into a double array.

If the input expression contains a symbolic variable, use VPA.

clear all
clc
% syms x
format longEng
syms y1 y2
h=4;
h1=2.88;
h2=0.26325;
alfa=0.5;
beta=0.5;
R1=-1; 
R3=-1;
q=0; 
kh=0.2;
kv=0.0;
psi=atan(kh/(1-kv));
gma=14; 
delta=20;
phi=30;
a=0.5;
b=1;


L=h+h1+h2;
psi=atan(kh/(1-kv));
da1=delta*(pi/180); pha1=phi*(pi/180)
l=pha1+da1;
m=pha1-psi;
n=psi+da1;

   B=(2*q*b*(b/(2*a+b)))/(gma*(h+h1)*(h+h1))
   alphacm=atan((sin(l)*sin(m)+(((sin(l)^2)*(sin(m)^2))-B*cos(n)*sin(m)*cos(l)+cos(m)*sin(l)*sin(m)*cos(l))^0.5)/((cos(m)*sin(l))-B*cos(n)))
   kg1=((tan(alphacm-phi)+tan(psi))*cos(alphacm-phi))/(tan(alphacm)*(cos(alphacm-phi-delta)))
   r1=(b/(h+h1))*(b/(2*a+b))*(tan(alphacm))
   kq1=r1*kg1
%    pgx=gma*(1-kv)*kg1*x
   pqx=(1-kv)*q*kq1
%

delm3=-0.5*(1-R3)*delta;
k3=(2*cos(phi-psi)^2)/(cos(phi-psi)^2*(1+R3)+cos(psi)*cos(delm3+psi)*(1-R3)*(1+sqrt((sin(phi+delm3)*sin(phi-psi))/cos(delm3+psi)))^2);


i=0;
for x=0:1:L
    i=i+1;
   %
%
        delm23=0.5*(3-1)*delta;
        k23=1+0.5*(3-1).*((cos(phi-psi).^2./(cos(psi).*(cos(delm23+psi).*...
            (-sqrt((sin(phi+delm23).*sin(phi-psi))./(cos(delm23+psi)))+1).^2)))-1);
       
        R2=3*(beta*(1-y1))^0.5;
        delm213=0.5*(R2-1)*delta;
        k213=1+0.5*(R2-1).*((cos(phi-psi).^2./(cos(psi).*(cos(delm213+psi).*...
            (-sqrt((sin(phi+delm213).*sin(phi-psi))./(cos(delm213+psi)))+1).^2)))-1);
    
        delm201=0.5*(1-R2).*delta;
        k201=(2*cos(phi-psi)^2)/(cos(phi-psi)^2*(1+R2)+cos(psi)*cos(delm201+psi)...
            *(1-R2)*(1+sqrt((sin(phi+delm201)*sin(phi-psi))/cos(delm201+psi)))^2);
        
         delm43=0.5*(3-1)*delta;
         k43=1+0.5*(3-1)*((cos(phi-psi)^2/(cos(psi)*(cos(delm43+psi)*(-sqrt((sin(phi+delm43)...
             *sin(phi-psi))/(cos(delm43+psi)))+1)^2)))-1)
        
         R4=3*(alfa*y2)^0.5;
         delm413=0.5*(R4-1)*delta;
         k413=1+0.5*(R4-1)*((cos(phi-psi)^2/(cos(psi)*(cos(delm413+psi)*(-sqrt((sin(phi+delm413)...
             *sin(phi-psi))/(cos(delm413+psi)))+1)^2)))-1)
         
         delm401=0.5*(1-R4).*delta;
         k401=(2*cos(phi-psi)^2)/(cos(phi-psi)^2*(1+R4)+cos(psi)*cos(delm401+psi)...
             *(1-R4)*(1+sqrt((sin(phi+delm401)*sin(phi-psi))/cos(delm401+psi)))^2);
    %
    H23=(k23*y1*cos(delm23));
    h23=gma*(x-h)^2*vpa(subs(H23,0,(1-(1/beta))))
    H213=(k213*cos(delm213)*y1);
    h213=gma*(x-h)^2*vpa(subs(H213,(1-(1/beta)),(1-(1/(9*beta)))))
    H201=(k201*y1*cos(delm201));
    h201=gma*(x-h)^2*vpa(subs(H201,(1-(1/(9*beta))),1));
    U=h23+h213+h201;
   %
    H123=(k23*y1*cos(delm23));
    h123=gma*h1^2*vpa(subs(H23,0,(1-(1/beta))));
    H1213=(k213*cos(delm213)*y1);
    h1213=gma*h1^2*vpa(subs(H213,(1-(1/beta)),(1-(1/(9*beta)))));
    H1201=(k201*y1*cos(delm201));
    h1201=gma*h1^2*vpa(subs(H201,(1-(1/(9*beta))),1));
    U1=h123+h1213+h1201;
   %
    H401=(k401*y2*cos(delm401));
    h401=gma*(x-h-h1)^2*vpa(subs(H401,0,(1/(9*alfa))));
    H413=(k413*y2*cos(delm413));
    h413=gma*(x-h-h1)^2*vpa(subs(H413,(1/(9*alfa)),(1/alfa)));
    H43=(k43*y2*cos(delm43));
    h43=gma*(x-h-h1)^2*vpa(subs(H43,(1/alfa),1));
    %
    Hb401=(k401*cos(delm401));
    hb401=(gma*(x+h)+q)*(x-h-h1)*vpa(subs(H401,0,(1/(9*alfa))));
    Hb413=(k413*cos(delm413));
    hb413=(gma*(x+h)+q)*(x-h-h1)*vpa(subs(H413,(1/(9*alfa)),(1/alfa)));
    Hb43=k43*cos(delm43);
    hb43=(gma*(x+h)+q)*(x-h-h1)*vpa(subs(H43,(1/alfa),1));
    %
    h4=h401+h413+h43+hb401+hb413+hb43;
    %
    h3=0.5*gma*k3*(x-h-h1)^2*cos(delta)+gma*x*k3*(x-h-h1)*cos(delta);
   %
    M23=(k23*y1*(1-y1)*cos(delta));
    m23=gma*(x-h)^3*vpa(subs(M23,0,(1-(1/beta))));
    M213=(k213*cos(delm213)*y1*(1-y1));
    m213=gma*(x-h)^3*vpa(subs(M213,(1-(1/beta)),(1-(1/(9*beta)))));
    M201=(k201*y1*cos(delm201)*(1-y1));
    m201=gma*(x-h)^3*vpa(subs(H201,(1-(1/(9*beta))),1));
    m2=m23+m213+m201;
% 
    M123=(k23*y1*(1-y1)*cos(delta));
    m123=gma*(x-h)^3*vpa(subs(M23,0,(1-(1/beta))));
    M1213=(k213*cos(delm213)*y1*(1-y1));
    m1213=gma*(x-h)^3*vpa(subs(M213,(1-(1/beta)),(1-(1/(9*beta)))));
    M1201=(k201*y1*cos(delm201)*(1-y1));
    m1201=gma*(x-h)^3*vpa(subs(H201,(1-(1/(9*beta))),1));
    m12=m123+m1213+m1201;
%
    M401=(k401*cos(delm401)*y2^2);
    m401=gma*(x-h-h1)^3*vpa(subs(M401,0,(1/(9*alfa))));
    M413=(k413*cos(delm413)*y2^2);
    m413=gma*(x-h-h1)^3*vpa(subs(M413,(1/(9*alfa)),(1/alfa)));
    M43=(k43*cos(delta)*y2^2);
    m43=gma*(x-h-h1)^3*vpa(subs(M43,(1/alfa),1));
    
    %
    Mb401=(k401*cos(delm401)*y2);
    mb401=(gma*(x+h)+q)*(x-h-h1)^2*vpa(subs(M401,0,(1/(9*alfa))));
    Mb413=(k413*cos(delm413)*y2);
    mb413=(gma*(x+h)+q)*(x-h-h1)^2*vpa(subs(M413,(1/(9*alfa)),(1/alfa)));
    Mb43=(k43*cos(delta)*y2);
    mb43=(gma*(x+h)+q)*(x-h-h1)^2*vpa(subs(M43,(1/alfa),1));
    %
    m4=m401+m413+m43+mb401+mb413+mb43;
    %
    m3=0.5*gma*k3*(x-h-h1)^3*cos(delta)*(2/3)+gma*x*k3*(x-h-h1)^2*cos(delta)*(1/2);
    %

    if(x<h)
        SF(i)=-0.5*(1-kv)*kg1*gma*x^2*cos(delta)-(1-kv)*q*kq1*x*cos(delta);
        SF1(i)=0;
        BM(i)=(1/6)*(1-kv)*kg1*gma*x^3*cos(delta)+(1-kv)*q*kq1*x*(x/2)*cos(delta);
        BM1(i)=0;
    elseif(x>=h && x<(h+h1))
        SF(i)=-0.5*kg1*(1-kv)*gma*x^2*cos(delta)-(1-kv)*q*kq1*x*cos(delta)+U;
        SF1(i)=0;
        BM(i)=(1/6)*kg1*(1-kv)*gma*x^3*cos(delta)+(1-kv)*q*kq1*x*(x/2)*cos(delta)-m2;
        BM1(i)=0;
    else
        SF(i)=-0.5*kg1*(1-kv)*gma*(h+h1)^2*cos(da1) + U1 -pqx*(h+h1)*cos(delta)-h4+h3;
        SF1(i)=0;
        BM(i)=0.5*kg1*(1-kv)*gma*(h+h1)^2*(((h+h1)/3)+(x-h-h1))*cos(da1)-m12+pqx*(h+h1)*(((h+h1)/2)+(x-h-h1))*cos(delta)-m4+m3;
        BM1(i)=0;  
    end
end
x=0:1:L;
f=max(BM)
 
subplot(2,1,1)
plot(x,SF,x,SF1)
grid on
xlabel('Length of the beam in m')
ylabel('Shear force  in KN/m')
title('Shear force diagram')
subplot(2,1,2)
plot(x,BM,x,BM1)
grid on
xlabel('Length of the beam in m')
ylabel('Bending Moment in KN-m/m')
title('Bending Moment diagram')

Answer 1

[Stephan]

During the first run of your loop U=0.0, which can be converted to double in line 144 and stored in SF(1). But in the second run with i=2 U is becomming a nonlinear equation depending on the symbolic variables y1 and y2. This expression can not be converted to double, because y1 and y2 dont have values:

U =        % Result of i = 1
 
0.0  
 
 
U =        % Result of i = 2
 
6.0365253422831727711781109416652*y1 + (0.050605831776812582566138498663122*y1*cos(30.0*(0.5 - 0.5*y1)^(1/2) - 10.0))/(0.0054220534046584909892291248567631*(0.5 - 0.5*y1)^(1/2) - 0.98058067569092015962081232865823*cos(30.0*(0.5 - 0.5*y1)^(1/2) - 10.197395559849880775082908712648)*(((0.99909591574838996219654063679627*sin(30.0*(0.5 - 0.5*y1)^(1/2) - 40.0))/cos(30.0*(0.5 - 0.5*y1)^(1/2) - 10.197395559849880775082908712648))^(1/2) + 1.0)^2*(3.0*(0.5 - 0.5*y1)^(1/2) - 1.0) + 0.0018073511348861636630763749522544) + 14*y1*cos(30.0*(0.38888888888888888888888888888889*y1 + 0.5)^(1/2) - 10.0)*(((0.0018431437409397227774065611925384*cos(30.0*(0.38888888888888888888888888888889*y1 + 0.5)^(1/2) - 9.8026044401501192249170912873524)^(7/9))/((-0.99909591574838996219654063679627*cos(30.0*(0.38888888888888888888888888888889*y1 + 0.5)^(1/2) - 9.8026044401501192249170912873524)^(7/9)*sin(30.0*(0.38888888888888888888888888888889*y1 + 0.5)^(1/2) + 20.0))^(1/2) + 0.77777777777777777777777777777778)^2 + 0.77777777777777777777777777777778)*(1.5*(0.38888888888888888888888888888889*y1 + 0.5)^(1/2) + 0.38888888888888888888888888888889) + 1.0)
 
The following error occurred converting from sym to double:
Unable to convert expression into double array.

Error in signal_check (line 145)
        SF(i)=-0.5*kg1*(1-kv)*gma*x^2*cos(delta)-(1-kv)*q*kq1*x*cos(delta)+U;

This is the point where things go wrong - you should think about how to fix this. I guess in line 75 where you calculate U will be a starting point for adapting your code:

U=h23+h213+h201;