# How to take integral of an equation?

2 views (last 30 days)
Ali Deniz on 9 Apr 2022
Commented: Torsten on 10 Apr 2022
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.

Torsten on 9 Apr 2022
Edited: Torsten on 9 Apr 2022
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);
##### 2 CommentsShowHide 1 older comment
Torsten on 10 Apr 2022
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)