How to take integral of an equation?

clc
clear all

hold on

R=0.56225;
delta_R=0.0015;
teta0=1.624*(pi/180);
teta1=90*(pi/180);
P1=0;
e=0.99;
P_8=zeros([],1);
teta_graph_8=linspace(0,90,8911);
W=zeros([],1);
i=1;

for P0=8
    
    for teta=-teta0:0.02*pi/80:teta0
        
        P_8(i)=8;
        i=i+1;
        
    end
    
    
    for teta=teta0:0.01*pi/180:teta1
        psi1=(e^2*(1-e^2)*cos(teta))/((1-e*cos(teta))^2)*(3-e^2-e/2*cos(teta)*(5-e^2))-e*(e^2+3)*log((1-e*cos(teta))/sin(teta))+(3*e^2+1)*log(tan(teta/2));
        psi2=(e^2*(1-e^2)*cos(teta0))/((1-e*cos(teta0))^2)*(3-e^2-(e)/2*cos(teta0)*(5-e^2))-e*(e^2+3)*log((1-e*cos(teta0))/sin(teta0))+(3*e^2+1)*log(tan(teta0/2));
        P_8(i)=(P1^2+(P0^2*(1-(P1/P0)^2)*(psi1/psi2)))^(1/2);
        syms teta
        W(i)=((pi*R^2)*(2*int((P1^2+(P0^2*(1-(P1/P0)^2)*(psi1/psi2)))^(1/2)*sin(teta)*cos(teta)+P0*(R0/R)^2-P1,teta,teta0,teta)));
        i=i+1;
        
        
    end
    
    
end

plot(teta_graph_8,P_8)
plot(-teta_graph_8,P_8)

W Formula.

I get the error "The following error occurred converting from sym to double:

Unable to convert expression into double array.

Error in Weight (line 32)

W(i)=((pi*R^2)*(2*int((P1^2+(P0^2*(1-(P1/P0)^2)*(psi1/psi2)))^(1/2),teta,teta0,teta))); " when I run this code. How can I take the formula?

Thank you.

Answer 1

[Torsten]

Your code is chaotic, so I don't know exactly what you are trying to do.

Make the best of this:

R=0.56225;
R0 = 0.1;
delta_R=0.0015;
theta0=1.624*(pi/180);
theta1=90*(pi/180);
P0=8;
P1=0;
e=0.99;
psi1 = @(theta) (e.^2.*(1-e.^2).*cos(theta))./((1-e.*cos(theta)).^2).*(3-e.^2-e./2.*cos(theta).*(5-e.^2))-e.*(e.^2+3).*log((1-e.*cos(theta))./sin(theta))+(3.*e.^2+1).*log(tan(theta./2));
psi2 = (e.^2.*(1-e.^2).*cos(theta0))./((1-e.*cos(theta0)).^2).*(3-e.^2-(e)./2.*cos(theta0).*(5-e.^2))-e.*(e.^2+3).*log((1-e.*cos(theta0))./sin(theta0))+(3.*e.^2+1).*log(tan(theta0./2));
P_8 = @(theta)(P1.^2+(P0.^2.*(1-(P1./P0).^2).*(psi1(theta)./psi2))).^(1/2);
fun = @(theta)P_8(theta).*sin(theta).*cos(theta)+P0*(R0./R).^2-P1;
W = 2*pi*R^2*integral(fun,theta0,theta1);

[Ali Deniz]

clc
clear all

hold on

R=0.56225;
R0=0.0015;
teta0=1.624*(pi/180);
teta1=90*(pi/180);
P1=0;
e=0.99;
P_8=zeros([],1);
teta_graph_8=linspace(0,90,8911);
W=zeros([],1);
i=1;

for P0=8
    
    for teta=-teta0:0.02*pi/80:teta0
        
        P_8(i)=8;
        i=i+1;
        
    end
    
    
    for teta=teta0:0.01*pi/180:teta1
        psi1=@(teta)(e^2*(1-e^2)*cos(teta))/((1-e*cos(teta))^2)*(3-e^2-e/2*cos(teta)*(5-e^2))-e*(e^2+3)*log((1-e*cos(teta))/sin(teta))+(3*e^2+1)*log(tan(teta/2));
        psi2=(e^2*(1-e^2)*cos(teta0))/((1-e*cos(teta0))^2)*(3-e^2-(e)/2*cos(teta0)*(5-e^2))-e*(e^2+3)*log((1-e*cos(teta0))/sin(teta0))+(3*e^2+1)*log(tan(teta0/2));
        P_8(i)=@(teta)(P1^2+(P0^2*(1-(P1/P0)^2)*(psi1/psi2)))^(1/2);
        
        fun=@(teta)((P1^2+(P0^2*(1-(P1/P0)^2)*(psi1/psi2)))^(1/2)*sin(teta)*cos(teta)+P0*(R0/R)^2-P1);
        W(i)=2*pi*R^2*integral(fun,teta0,teta1);
        i=i+1;
        
        
    end
    
    
end

plot(teta_graph_8,P_8)
plot(-teta_graph_8,P_8)
plot(P_8,W);

I did this but now the error "Conversion to double from function_handle is not possible.

Error in Weight (line 30)

P_8(i)=@(teta)(P1^2+(P0^2*(1-(P1/P0)^2)*(psi1/psi2)))^(1/2);" . Thank you.

[Torsten]

P0 = 8;
P1 = 0;
e = 0.99;
theta0 = 1.624*(pi/180);
theta1 = 90*(pi/180);
R0 = 0.0015;
R = 0.56225;
psi1 = @(theta) (e.^2.*(1-e.^2).*cos(theta))./((1-e.*cos(theta)).^2).*(3-e.^2-e./2.*cos(theta).*(5-e.^2))-e.*(e.^2+3).*log((1-e.*cos(theta))./sin(theta))+(3.*e.^2+1).*log(tan(theta./2));
psi2 = (e.^2.*(1-e.^2).*cos(theta0))./((1-e.*cos(theta0)).^2).*(3-e.^2-(e)./2.*cos(theta0).*(5-e.^2))-e.*(e.^2+3).*log((1-e.*cos(theta0))./sin(theta0))+(3.*e.^2+1).*log(tan(theta0./2));
P_8 = @(theta)(P1.^2+(P0.^2.*(1-(P1./P0).^2).*(psi1(theta)./psi2))).^(1/2);
fun = @(theta)P_8(theta).*sin(theta).*cos(theta);
theta = theta0:0.01*pi/180:theta1;
W = zeros(numel(theta),1);
for i = 1:numel(theta)-1
    theta_low = theta(i);
    theta_high = theta(i+1);
    value_integral = integral(fun,theta_low,theta_high);
    W(i+1) = W(i) + value_integral;
end
W = pi*R^2*(2*W + P0*(R0/R)^2 - P1);
plot(theta,W)