Problems using diff() to compute function handle derivative

Hi everybody,

First of all, sorry for the large code-insert. But as explained below, I wasn't able to reproduce the error within a simple example.

I try to compute the derivative of the function handles "B_Teeth" and "B_Yoke" with respect to the only variable t (see last few rows of the code). As I need to transform the output of this operation further, I need to have it as a function of t itself.

As I read here, this is in principle feasable but although I've been trying for a while it does not work. What I'm wondering about are two things:

1) The error code states that the input argument of diff() is sym. But when I doublecheck this via class(), it says that "B_Teeth" is a function handle (as it should be).

2) I tried to simplify the problem and reproduce the error code (see second code) but in this case, everything works well.

I highly appreciate your help!

%% Parameter
p=4;N=32;a_p=0.742;b=0.886;B_r=1.25;mu_r=1.024;Rs=0.091;Rr=0.090;t_b=0.001;w_m=0.020;w_r=0.008;mu_0=4*pi*10^(-7);   
l_m=0.006;l_bar1=0.005;l_bar2=0.0059;w_bar=0.005;w_b=w_bar+2*t_b; Q=48;L=135/149*Rr;Wtau_p=5/6;m=3;q=Q/(2*p*m);J=10;Omega=100;a_PM=140/360*2*pi;
tau_p=pi*Rs/p;tau_s=pi*2*Rs/Q;b_StZa=0.004;a_s=pi/(m*q);mu_b=2.3696e-5;tau_2=pi*(Rs+0.004)/p; B_b=1.6178;
phi=pi; I=100; k_v=[-1;1;-1;1;-1;1]; v=[1;5;7;11;13;17];
%% Reaktances
P1=mu_0./(p.*(Rs-Rr)).*((Rs+Rr)./2).*L; P2=P1.*a_p.*pi; P3=2.*mu_0.*mu_r.*w_m.*L./l_m; P4=2*mu_0*(l_bar1+l_bar2)*L/w_bar; P5=2*mu_b*t_b*L/w_b;
P6=mu_0/(p*(Rs-Rr))*(Rs+Rr)/2*L*pi*(b-a_p)/2;

syms nu theta j t k_vv
k_wv=sin(Wtau_p.*nu.*pi./2).*(sin(nu.*pi./(2.*m))./(q.*sin(nu.*pi./(2.*m.*q))));
F_rdModHiVarUd1=4./pi.*symsum(sin((2*j+1).*a_p.*pi./2)./(2*j+1).*cos((2*j+1).*(theta-Omega.*t)),j,0,J);
F_rdModHiVarUd2=4./pi.*symsum((sin((2*j+1).*b.*pi./2)-sin((2*j+1).*a_p.*pi./2))./(2*j+1).*cos((2.*j+1).*(theta-Omega.*t)),j,0,J); 
Fr_qModHiVar=4./pi.*symsum((cos((2*j+1).*a_p.*pi./2)-cos((2*j+1).*b.*pi./2))./(2*j+1).*sin((2.*j+1).*(theta-Omega.*t)),j,0,J);
B_airgapPM=B_r.*4./pi.*symsum(sin((2*j+1).*a_p.*pi./2)./(2.*j+1).*cos((2.*j+1).*(theta-Omega.*t)),j,0,J);
F=3.*N.*I.*k_wv./(nu.*p.*pi); 
U_d1Mod_HiVar=2.*P1.*F.*sin(nu.*a_p.*pi./2)./(nu.*(P2+P3+P4+P5)).*cos(phi).*cos((k_vv+nu).*Omega.*t);
U_d2Mod_HiVar=2.*P1.*F.*sin(nu.*pi./4.*(b-a_p)).*cos(phi)./(nu.*(0.5.*P5+P6)).*cos((nu+k_vv).*Omega.*t+nu.*(b+a_p).*pi./4);
U_qMod_HiVar=-2.*k_vv.*P1.*F.*sin(phi).*sin(nu.*pi./4.*(b-a_p)).*sin(Omega.*t.*(k_vv+nu).*nu+pi./4.*(a_p+b))./(nu.*(2.*P5+P6));
for i=1:length(v)
    U_d1Mod(i)=subs(U_d1Mod_HiVar,[nu,k_vv],[v(i),k_v(i)]);
    U_d2Mod(i)=subs(U_d2Mod_HiVar,[nu,k_vv],[v(i),k_v(i)]);
    U_qMod(i)=subs(U_qMod_HiVar,[nu,k_vv],[v(i),k_v(i)]);
end
U_d1Mod=sum(U_d1Mod,'all');
U_d2Mod=sum(U_d2Mod,'all');
U_qMod=sum(U_qMod,'all');
Fs_dMod=F.*cos(nu.*theta+k_vv.*Omega.*t).*cos(phi); 
Fs_qMod=F.*k_vv.*sin(-k_vv*Omega.*t-nu.*theta).*sin(phi);
Fr_dMod=U_d1Mod.*F_rdModHiVarUd1+U_d2Mod.*F_rdModHiVarUd2;
Fr_qMod=U_qMod.*Fr_qModHiVar;
for i=1:length(v)
Fs_DMod(i)=subs(Fs_dMod,[nu,k_vv],[v(i),k_v(i)]);
Fs_QMod(i)=subs(Fs_qMod,[nu,k_vv],[v(i),k_v(i)]);
end
Fs_DMod=(sum(Fs_DMod,'all'));
Fs_QMod=(sum(Fs_QMod,'all'));
B_dMod=(Fs_DMod-Fr_dMod).*mu_0./(Rs-Rr);
B_qMod=(Fs_QMod-Fr_qMod).*mu_0./(Rs-Rr);

%% Airgap Flux density in yoke and teeth 
B_airgapMod=matlabFunction(B_dMod+B_qMod+B_airgapPM);
syms t 
B_Teeth=@(t)tau_s./(b_StZa.*a_s).*integral(@(theta)B_airgapMod(t,theta),-a_s/2,a_s/2);
B_Yoke=@(t)tau_p./(3.*tau_2.*pi).*sqrt(2./(3.*q)).*integral(@(theta)B_airgapMod(t,theta),-pi/2,pi/2);

DeltaB_Yoke=matlabFunction(diff(B_Yoke(t)));
DeltaB_Teeth=matlabFunction(diff(B_Teeth(t)));

The good-working simple trial is:

syms x
Fun=@(x)cos(x);
DiffFun=diff(Fun(x));
FunDiffFun=matlabFunction(DiffFun);
FunDiffFun(2)

Answer 1

[Walter Roberson]

You do not have any choice. You need a symbolic expression of B_teeth and you cannot get one by doing numeric integration. You will have to do symbolic integration.

vpaintegral() will only evaluate to a placeholder if there are any free variables not being integrated over. However, diff() is able to reason about those placeholders, and so there is a possibility that you might be able to get somewhere. I would not count on it though.

Note: if you do use int() on a summation, in some cases it is faster to take children() of the summation to get the individual terms, int() those, and add the results together. This will keep int() from trying to pattern match each combination of terms to find something in its library that it has a rule for, which can be a waste of time if you have reasons to believe that there is no such useful combination.

[Richard]

I will both try out the possibility mentioned by Star Strider (to use vpaintegral()) and yours (to use children()).

I'll let you know what worked or not!

Anyway, sincere thanks to both for advising me!

[Richard]

As promised I would like to document briefly how I have solved the problem. I have stuck to the symbolic integration as I was advised. I could avoid solving the integral in every iteration loop by defining also those values as symbolic variables, which change between two iterations, but are constant within one iteration step. (E.g. speed omega is such a variable, while time t could take different values even within one iteration step). This integral with only general relations, I can then derive symbolically and process further. During the iteration I then only substitute all variables and solve the prefabricated expressions for the integral and the derivation. The disadvantage of this method is, that already the insertion of the numerical values into the symbolic expressions and the calculation via eval(subs(.)) takes several seconds with a high number of Fourier coefficients. However, this duration is still acceptable for my purposes, so I have stuck to it.

Since this solution works for me and I have time pressure, I could not yet test the proposed alternatives via children().

Finally, I woud like to thank you both very much for the support and the invested time! The help of both of you has been extremly helpful, so best would be to accept both your answers. This is not possible, so I accept the last answer without any less appreciation for the work of Star Strider!