How to vary a value inside an integral which is not getting integrated itself?

I am trying to find the triple integral of a function with limits defined for 'phi', 'z' and 'r'. On finding the value of the integral at any particular value say h=20 (h IS NOT GETTING INTEGRATED), I get a correct answer which matches with my experimental values. However, as soon as I vary (using linspace or colon) the value of h (h IS NOT GETTING INTEGRATED) I get an error saying "Matrix dimensions must agree". I have attached the code below for reference. Thank you

clc,clear

%magnetic field for an off axis point due to a thick solenoid in z-x plane

l=20; %length of the solenoid

r1 = 5; %inner radius

r2 = 10; %outer radius

u0 = 1.26*(10^-3); %permeability of free space

i = 3; %current in amps

N = 800; %number of turns

C=(u0*i*N)/(2*pi*10*l);

%location of x,y,z cooridnates of P

%rho = linspace(0,22,100);

rho = 0;

h=20:1:35;

Binz = @(r,z,phi)(r.*(r-(rho.*cos(phi))))./(((rho.^2) + (r.^2) + ((h-z).^2) - (2.*r.*rho.*cos(phi))).^(3/2));

Bz = C*integral3(Binz,r1,r2,0,20,0,2*pi); %magnetic field in z direction

plot (Bz,h);

Answer 1

[Star Strider]

The integral2 and integral3 functions do not allow the 'ArrayValued' name-value pair, so it is necessary to ‘nest’ integral calls if you want to do that.

Try this:

%magnetic field for an off axis point due to a thick solenoid in z-x plane
l=20;                                             %length of the solenoid
r1 = 5;                                           %inner radius 
r2 = 10;                                         %outer radius
u0 = 1.26*(10^-3);                        %permeability of free space
i = 3;                                             %current in amps
N = 800;                                        %number of turns
C=(u0*i*N)/(2*pi*10*l);
%location of x,y,z cooridnates of P
%rho = linspace(0,22,100);
rho = 0;
h=20:1:35;
Binz = @(r,z,phi)(r.*(r-(rho.*cos(phi))))./(((rho.^2) + (r.^2) + ((h-z).^2) - (2.*r.*rho.*cos(phi))).^(3/2));
% Bz = C*integral3(Binz,r1,r2,0,20,0,2*pi);    %magnetic field in z direction
Bz = integral(@(phi) integral(@(z) integral(@(r) Binz(r,z,phi), r1,r2, 'ArrayValued',1), 0,20, 'ArrayValued',1),0,2*pi, 'ArrayValued',1)
plot (Bz,h);

producing:

Bz =
   29.3733   25.2358   21.3369   17.8589   14.8792   12.3925   10.3485    8.6806    7.3222    6.2144    5.3077    4.5619    3.9451    3.4317    3.0019    2.6398

and this plot:

Please check that to be certain it appears reasonable.

Just to be sure, I checked this approach with the first example in the integral3 documentation:

fun = @(x,y,z) y.*sin(x)+z.*cos(x)
q = integral(@(z) integral(@(y) integral(@(x) fun(x,y,z), 0, pi, 'ArrayValued',1), 0,1, 'ArrayValued',1), -1,1, 'ArrayValued',1)

with:

q =
    2.0000

being the same as the result in the integral3 documentation.

[Azeem Singh Kahlon]

Thanks a lot for the answer the values perfectly match. Thanks for the help.

[Star Strider]

As always, my pleasure!