Problema con syms y el calculo de integrales
Publicado por raimundo (1 intervención) el 08/06/2011 06:08:58
Estimados usuarios,
Soy nuevo en matlab, y estoy aprendiendo a usarlo debido a una tarea que tengo que hacer para la universidad. Se trata de hacer un diagrama de momento curvatura para una seccion doblemente armada de hormigon. Fuera de esto, queria saber si porcasualidad alguien me puede ayudar, ya que dentro del archivo tipo .m tengo unas integrales metidas que las queria calcular de la forma syms x y despues int(funcion,a,b) pero el programa me tira el siguiente error:
??? Undefined function or method 'le' for input arguments of type 'sym'.
Error in ==> DMC at 134
if F<=2 & F>=-2
adjunto el codigo escrito y agradeceria demasiado si alguien me pudiese ayudar.
%Momento Curvatura De una Seccion doblemente armada
syms y
disp('DIAGRAMA DE MOMENTO CURVATURA DE UNA SECCION RECTANGULAR DOBLEMENTE ARMADA')
disp('')
disp(' El siguiente programa supone siempre el uso de acero A630-420H y hormigon H30')
b=input('Introduzca el ancho de la sección en centimetros ');
h=input('La altura de la seccion en cm ');
a1=input('El area de acero As inferior en cm2 ');
a2=input('El area de acero As prima (superior) en cm2 ');
r1=input('El recubrimiento inferior en cm ');
r2=input('El recubrimiento prima (superior) en cm ');
Es=2.1*10^6; %Mod elasticidad acero en kgf/Cm2
%Primer punto del grafico corresponde al agrietamiento por traccion del
%hormigon. Se obtiene con resistencia de materiales (navier y usando el fc de traccion) y
%suponiendo que c=h/2 (aproximacion)
I=(b*h^3)/12; %momento de inercia
E=4700*sqrt(25); %modulo de elasticidad del hormigon en MPa , segun ACI318
ft=0.62*sqrt(25); %resistencia a traccion en MPa, segun ACI318
%Entonces usando la ecuacion de navier se obtiene el momento de resistencia
%a tracción y la curvatura en ese momento con la expresion fi = M/EI
%Luego
m1=ft*10*I/(h/2); %El factor 10 es para pasar de MPa a kgf/cm2. En la siguiente linea del codigo tambien se usa el factor 10.
f1=m1/(E*10*I);
fi(1)=0; %El vector fi y el vector mo seran usados para graficar el resultado final que se busca. El diagrama M vs fi (curvatura);
mo(1)=0;
%El grafico parte del punto (0,0)
fi(2)=f1;
mo(2)=m1; %Segundo punto del grafico, correspendiente a la falla a tracción del hormigon
%Se tiene asi el primer punto de la curva
%Luego hay que iterar c, para que se cumpla la compatibilidad de
%deformaciones y el equilibrio para el siguiente tramo del grafico.
%Se desprecia el aporte del hormigon como fuerza en traccion, y se da un
%valor inicial para C (distancia desde el borde sup hasta el eje neutro).
%Tambien se da un valor inicial para ec muy pequeño que este a priori
%cercano al ec q se cumple cuando falla a traccion) ec es la deformacion
%unitaria del hormigon en su borde superior.
j=2;
for ec=0.0003:0.00015:0.003 %Hasta donde falla el hormigon en la parte superior
j=1+j; %Contador para generar un vector que guarde los valores a graficar
%De esta manera se parte con j=3, es decir el tercer punto del grafico y
%asi sucesivamente
c=0.5; %Valor inicial para c , de donde se empieza a iterar
F=3;
while c<h && abs(F)>2
c=c+0.5; %aumento c a razon de 0,5 cada vez
es1=ec*(h-c-r1)/c;
es2=ec*(c-r2)/c; %compatibilidad de deformaciones, es1 es la def del acero inferior y es2 del superior.
%Se utiliza el modelo Park para la curva esf. deformacion del hormigon
%, considerando los limites de integracion y una relacion entre ec
% ,la cantidad c y la altura y medida desde el eje neutro hasta el borde
% superior, que se obtiene con geometria. esto es: ey=ec*y/c
if es1>0 && es2>0 %Para usar solo valores fisicamente validos
if es1<0.002 & es2<0.002 & ec<=0.002
t1=Es*es1*a1; %t1 es la fuerza de tracción hecha por el acero inferior
c2=Es*es2*a2; %c2 es la fuerza a compresion hecha por el acero superior
fc=b*250*(ec*c/0.002-(c/3)*(ec/0.002)^2); %Fza a compresion del hormigon
end
if es1>0.002 & es2<0.002 & ec<=0.002
t1=4200*a1;
c2=Es*es2*a2;
fc=b*250*(ec*c/0.002-(c/3)*(ec/0.002)^2);
end
if es1>0.002 & es2<0.002 & ec>0.002
t1=4200*a1;
c2=Es*es2*a2;
fc=b*250*(((2/3)*0.002*c/ec)+int(1-100*(ec*y/c-0.002),0.002*c/ec,c));
end
if es1>0.002 & es2>0.002 & ec>0.002
t1=4200*a1;
c2=4200*a2;
fc=b*250*(((2/3)*0.002*c/ec)+int(1-100*(ec*y/c-0.002),0.002*c/ec,c));
end
if es1<0.002 & es2<0.002 & ec>0.002
t1=Es*es1*a1;
c2=Es*es2*a2;
fc=b*250*(((2/3)*0.002*c/ec)+int(1-100*(ec*y/c-0.002),0.002*c/ec,c));
end
if es1>0.002 & es2>0.002 & ec<=0.002
t1=4200*a1;
c2=4200*a2;
fc=b*250*(ec*c/0.002-(c/3)*(ec/0.002)^2);
end
F=t1-c2-fc;
if F<=2 & F>=-2
fi(j)=ec/c;
if ec<=0.002
mo(j)=t1*(h-r1-c)+c2*(c-r2)+ b*250*int(2*y*(ec.*y/0.002*c)-y*(ec.*y/0.002*c)^2,0,c); %torque resp. al eje neutro
end
if ec>0.002
mo(j)=t1*(h-r1-c)+c2*(c-r2)+b*250*(int(2*y*(ec*y/0.002*c)-y*(ec*y/0.002*c)^2,0,(0.002*c/ec))+int(y-100*y*(ec*y/c-0.002),(0.002*c/ec),c));
end
end
end %Corresponde a la condicion de deformaciones positivas
end
end
plot(fi,mo)
Soy nuevo en matlab, y estoy aprendiendo a usarlo debido a una tarea que tengo que hacer para la universidad. Se trata de hacer un diagrama de momento curvatura para una seccion doblemente armada de hormigon. Fuera de esto, queria saber si porcasualidad alguien me puede ayudar, ya que dentro del archivo tipo .m tengo unas integrales metidas que las queria calcular de la forma syms x y despues int(funcion,a,b) pero el programa me tira el siguiente error:
??? Undefined function or method 'le' for input arguments of type 'sym'.
Error in ==> DMC at 134
if F<=2 & F>=-2
adjunto el codigo escrito y agradeceria demasiado si alguien me pudiese ayudar.
%Momento Curvatura De una Seccion doblemente armada
syms y
disp('DIAGRAMA DE MOMENTO CURVATURA DE UNA SECCION RECTANGULAR DOBLEMENTE ARMADA')
disp('')
disp(' El siguiente programa supone siempre el uso de acero A630-420H y hormigon H30')
b=input('Introduzca el ancho de la sección en centimetros ');
h=input('La altura de la seccion en cm ');
a1=input('El area de acero As inferior en cm2 ');
a2=input('El area de acero As prima (superior) en cm2 ');
r1=input('El recubrimiento inferior en cm ');
r2=input('El recubrimiento prima (superior) en cm ');
Es=2.1*10^6; %Mod elasticidad acero en kgf/Cm2
%Primer punto del grafico corresponde al agrietamiento por traccion del
%hormigon. Se obtiene con resistencia de materiales (navier y usando el fc de traccion) y
%suponiendo que c=h/2 (aproximacion)
I=(b*h^3)/12; %momento de inercia
E=4700*sqrt(25); %modulo de elasticidad del hormigon en MPa , segun ACI318
ft=0.62*sqrt(25); %resistencia a traccion en MPa, segun ACI318
%Entonces usando la ecuacion de navier se obtiene el momento de resistencia
%a tracción y la curvatura en ese momento con la expresion fi = M/EI
%Luego
m1=ft*10*I/(h/2); %El factor 10 es para pasar de MPa a kgf/cm2. En la siguiente linea del codigo tambien se usa el factor 10.
f1=m1/(E*10*I);
fi(1)=0; %El vector fi y el vector mo seran usados para graficar el resultado final que se busca. El diagrama M vs fi (curvatura);
mo(1)=0;
%El grafico parte del punto (0,0)
fi(2)=f1;
mo(2)=m1; %Segundo punto del grafico, correspendiente a la falla a tracción del hormigon
%Se tiene asi el primer punto de la curva
%Luego hay que iterar c, para que se cumpla la compatibilidad de
%deformaciones y el equilibrio para el siguiente tramo del grafico.
%Se desprecia el aporte del hormigon como fuerza en traccion, y se da un
%valor inicial para C (distancia desde el borde sup hasta el eje neutro).
%Tambien se da un valor inicial para ec muy pequeño que este a priori
%cercano al ec q se cumple cuando falla a traccion) ec es la deformacion
%unitaria del hormigon en su borde superior.
j=2;
for ec=0.0003:0.00015:0.003 %Hasta donde falla el hormigon en la parte superior
j=1+j; %Contador para generar un vector que guarde los valores a graficar
%De esta manera se parte con j=3, es decir el tercer punto del grafico y
%asi sucesivamente
c=0.5; %Valor inicial para c , de donde se empieza a iterar
F=3;
while c<h && abs(F)>2
c=c+0.5; %aumento c a razon de 0,5 cada vez
es1=ec*(h-c-r1)/c;
es2=ec*(c-r2)/c; %compatibilidad de deformaciones, es1 es la def del acero inferior y es2 del superior.
%Se utiliza el modelo Park para la curva esf. deformacion del hormigon
%, considerando los limites de integracion y una relacion entre ec
% ,la cantidad c y la altura y medida desde el eje neutro hasta el borde
% superior, que se obtiene con geometria. esto es: ey=ec*y/c
if es1>0 && es2>0 %Para usar solo valores fisicamente validos
if es1<0.002 & es2<0.002 & ec<=0.002
t1=Es*es1*a1; %t1 es la fuerza de tracción hecha por el acero inferior
c2=Es*es2*a2; %c2 es la fuerza a compresion hecha por el acero superior
fc=b*250*(ec*c/0.002-(c/3)*(ec/0.002)^2); %Fza a compresion del hormigon
end
if es1>0.002 & es2<0.002 & ec<=0.002
t1=4200*a1;
c2=Es*es2*a2;
fc=b*250*(ec*c/0.002-(c/3)*(ec/0.002)^2);
end
if es1>0.002 & es2<0.002 & ec>0.002
t1=4200*a1;
c2=Es*es2*a2;
fc=b*250*(((2/3)*0.002*c/ec)+int(1-100*(ec*y/c-0.002),0.002*c/ec,c));
end
if es1>0.002 & es2>0.002 & ec>0.002
t1=4200*a1;
c2=4200*a2;
fc=b*250*(((2/3)*0.002*c/ec)+int(1-100*(ec*y/c-0.002),0.002*c/ec,c));
end
if es1<0.002 & es2<0.002 & ec>0.002
t1=Es*es1*a1;
c2=Es*es2*a2;
fc=b*250*(((2/3)*0.002*c/ec)+int(1-100*(ec*y/c-0.002),0.002*c/ec,c));
end
if es1>0.002 & es2>0.002 & ec<=0.002
t1=4200*a1;
c2=4200*a2;
fc=b*250*(ec*c/0.002-(c/3)*(ec/0.002)^2);
end
F=t1-c2-fc;
if F<=2 & F>=-2
fi(j)=ec/c;
if ec<=0.002
mo(j)=t1*(h-r1-c)+c2*(c-r2)+ b*250*int(2*y*(ec.*y/0.002*c)-y*(ec.*y/0.002*c)^2,0,c); %torque resp. al eje neutro
end
if ec>0.002
mo(j)=t1*(h-r1-c)+c2*(c-r2)+b*250*(int(2*y*(ec*y/0.002*c)-y*(ec*y/0.002*c)^2,0,(0.002*c/ec))+int(y-100*y*(ec*y/c-0.002),(0.002*c/ec),c));
end
end
end %Corresponde a la condicion de deformaciones positivas
end
end
plot(fi,mo)
Valora esta pregunta
0
