How to integrate an array properly in matlab

30 次查看(过去 30 天)
Hello, I am trying to finish a m-file to find the inductance in 2 coils. I have finished the program to the point of integration. I believe the problem is the integration of an array. I have tried various methods such as int, trapz, and quad but all of these seem to be returning an error. I am not sure as to whether I am implementing the commands wrong or rather I have a bad equation. Here is my code
% This program finds mutual inductance for TET coil system %
% Using Neumann's definition %
% M = sqrt(ap*as)*(1/(2*pi))*(int(((cos(theta)-(d/as)*((cos(psi)*cos(phi))-(sin(psi)*sin(phi)*cos(theta))))/(R^(3/2)))*f,phi,0,2*pi))
%
% Increment angles
%psi = [0:1:90];
%theta = 0:1:90;
%phi = 0:1:90;
psi = input('Enter psi value in degrees \n')
theta = input('Enter theta value in degrees \n')
phi = input('Enter phi value in degrees \n')
% Arguments
%dim = 0:2*pi
h = input('Enter h value in mm \n')
ap = input('Enter ap value in mm \n')
as = sqrt((ap^2)-(h^2))
delta = h/ap
alpha = as/ap
d = h.*tan(phi)
Ra = (1-(cos(phi).*cos(phi).*sin(theta).*sin(theta)))
Rb = ((2)*(d/as)).*((sin(psi).*sin(phi))-(cos(psi).*cos(phi).*cos(theta)))
Rc = ((d.^2)/(as.^2))
R = sqrt(Ra+Rb+Rc)
z = delta-(alpha*sin(theta)*cos(phi))
%kprime_2 = (((1-(alpha.*R)).^2)+z.^2)/(((1+(alpha.*R)).^2+z.^2))
kprime_2a = ((1-(alpha.*R)).^2)+z.^2
kprime_2b = ((1+(alpha.*R)).^2)+z.^2
kprime_2 = kprime_2a./kprime_2b
f = -0.011*(log(kprime_2))-0.0021
integrand = ((cos(theta)-(d/as).*((cos(psi).*cos(phi))-
(sin(psi).*sin(phi).*cos(theta))))/(R.^(3/2))).*f
%integrand1 = double(int(((cos(theta)-(d/as).*((cos(psi).*cos(phi))-(sin(psi).*sin(phi).*cos(theta))))./(R.^(3/2))).*f,phi,0,2.*pi))
%integrand = double(integrand);
stuff = trapz(phi,integrand)
M = sqrt(ap.*as).*(1/(2*pi)).*stuff
I am setting psi and phi to 0 and setting theta to 0:10:90. H is usually 3 and ap is usuall 6. This gives me different error messages for each method of integration I use.
Any help would be appreciated. Thanks.
  5 个评论
Andrew Matthews
Andrew Matthews 2013-6-7
编辑:Andrew Matthews 2013-6-7
Whenever I use 1 my integral becomes 0. It was not supposed to be 0. Would you know the reasoning? Thank you.
Andrew Matthews
Andrew Matthews 2013-6-7
I just realized, it is because my program is filling in the angles before integration. Because of this I am getting 0.

请先登录,再进行评论。

采纳的回答

Andrew Newell
Andrew Newell 2013-6-6
编辑:Andrew Newell 2013-6-6
Assuming you are trying to integrate over theta, a good approach is to create a function for your integrand like this:
function y = inductanceIntegrand(theta,psi,phi,h,ap)
as = sqrt((ap^2)-(h^2));
delta = h/ap;
alpha = as/ap;
d = h.*tan(phi);
Ra = (1-(cos(phi).*cos(phi).*sin(theta).*sin(theta)));
Rb = ((2)*(d/as)).*((sin(psi).*sin(phi))-(cos(psi).*cos(phi).*cos(theta)));
Rc = ((d.^2)/(as.^2));
R = sqrt(Ra+Rb+Rc);
z = delta-(alpha*sin(theta).*cos(phi));
kprime_2a = ((1-(alpha.*R)).^2)+z.^2;
kprime_2b = ((1+(alpha.*R)).^2)+z.^2;
kprime_2 = kprime_2a./kprime_2b;
f = -0.011*(log(kprime_2))-0.0021;
y = sqrt(ap.*as).*(1/(2*pi)).*((cos(theta)-(d/as).*((cos(psi).*cos(phi))-(sin(psi).*sin(phi).*cos(theta))))/(R.^(3/2))).*f;
Then you can create an anonymous function that is a parameter of theta only, and integrate:
psi = 0;
phi = 0;
h = 3;
ap = 6;
f = @(theta) inductanceIntegrand(theta,psi,phi,h,ap);
M = quadl(f,0,2*pi);
disp(M)
-4.5074e-04
If you're integrating over phi, just change the argument of f to phi and provide a value for theta.
  3 个评论
Andrew Matthews
Andrew Matthews 2013-6-7
编辑:Andrew Matthews 2013-6-7
When I run the code the way you have it written it says 'requires more input arguments to run'. And when I have variables for input when I run the second m file it says the output and input vectors must be the same length. What do I do from here? Thank you.
Andrew Matthews
Andrew Matthews 2013-6-7
I just realize my error. I was trying to run to man variables at a time. Thank you all, you have been a great help.

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by