Symbolic derivative of an arrayfun

2 次查看（过去 30 天）
Maurilio Matracia
Maurilio Matracia 2021-4-18
Hello everyone,
I am trying to compute the derivative of the piecewise function F_Ma, which I called f_Ma.
infty = 1e9; zero = 1e-9 ; Pi = 3.1415 ;
%% PARAMETERS
type =1;
p = [1.585 1.585 10] ;
H=100 ; NA = 10 ;
urbLev = 3 ;
alpha = [2, 3, 3] ;
tau = 0.31623 ; sigma_n2 = 1e-12 ;
m = [2 1 1] ;
Rd = 1e3 ;
Ru = [0 0.3 0.7]' * Rd ;
eta = [0.9772 0.007943 1 ; 0.7943 0.01 1 ; 0.6918 0.005 1 ; 0.5888 0.0004 1] ; eta(:,3) = eta(:,1) ;
Sa = [4.88 9.6117 12.0810 27.304]' ; Sb = [0.429 0.1581 0.1139 0.0797]' ;
nsteps = 1 ;
Z = linspace(zero, max(Ru)+Rd, nsteps) ;
%% CDF ANALYSIS
for u = 1:length(Ru)
IND{u,1} = @(z) z<= Rd-Ru(u) ;
IND{u,2} = @(z) (Rd-Ru(u) < z) & (z < Rd+Ru(u)) ;
IND{u,3} = @(z) z >= Rd+Ru(u) ;
Betai{u} = @(z) real( (z>0).*acos( (z.^2 + Ru(u).^2 - Rd.^2) ./ ( (z>0).*2.*z.*Ru(u) + (z==0)) ) ) ;
dBetai_dz{u} = @(z) (z>0 ).*(Ru(u).^2 - z.^2 - Rd.^2) ./ ...
( (z>0) .* (z .* sqrt(4*z.^2 .* Ru(u).^2 - (z.^2 + Ru(u).^2 - Rd.^2).^2)) + (z==0) ) ;
Betai{u} = @(z) real( acos( (z.^2 + Ru(u).^2 - Rd.^2) ./ ( 2.*z.*Ru(u) ) ) ) ;
dBetai_dz{u} = @(z) (Ru(u).^2 - z.^2 - Rd.^2) ./ ...
( z .* sqrt(4*z.^2 .*Ru(u).^2 - (z.^2 + Ru(u).^2 - Rd.^2).^2) ) ;
zX{u} = @(beta) Ru(u).*cos(beta) + sqrt(Rd.^2 - Ru(u).^2.*sin(beta).^2) ;
Sigma{1} = @(z) pi*z.^2 ;
Sigma{2} = @(z) Betai{u}(z).*z.^2 + 1/2*integral(@(beta) zX{u}(beta).^2, Betai{u}(z),2*pi-Betai{u}(z),'ArrayValued',1) ;
Sigma{3} = @(z) pi*Rd^2 ;
dSigma{1} = @(z) 2*pi*z ;
dSigma{2} = @(z) 2*z.*Betai{u}(z) + z.^2 .* dBetai_dz{u}(z) - zX{u}(Betai{u}(z)).^2 .* dBetai_dz{u}(z) ;
dSigma{3} = @(z) 0 ;
P_Ma{1} = @(z) 1./(1+Sa(urbLev)*exp(-Sb(urbLev)*(180/pi*atan(H./z)-Sa(urbLev)))) ;
P_Ma{2} = @(z) 1-P_Ma{1}(z) ;
for M = 1:2
for i = 1:3
fDelta_i = @(d,z) 1./Sigma{i}(z) .* (dSigma{1}(d).*IND{u,1}(d)+dSigma{2}(d).*IND{u,2}(d)+dSigma{3}(d).*IND{u,3}(d)); % Only because we are going to integrate over d from 0 to z
Ups_i = @(z_) arrayfun (@(z) integral(@(d) fDelta_i(d,z) .* (1-P_Ma{M}(d)), zero,z,'ArrayValued',1), z_) ;
dUps_i = @(z_) arrayfun(@(z) fDelta_i(z,z).*(1-P_Ma{M}(z)) - dSigma{i}(z)./(Sigma{i}(z).^2) .* ...
integral(@(d) dSigma{1}(d).*IND{u,1}(d)+dSigma{2}(d).*IND{u,2}(d)+dSigma{3}(d).*IND{u,3}(d), 0,z,'ArrayValued',1), z_) ;
ps = @(z) Sigma{i}(z) / (pi*Rd^2) ;
dps = @(z) dSigma{i}(z) / (pi*Rd^2) ;
P_k = @(z,k) nchoosek(NA,k) .* ps(z).^k .* (1-ps(z)).^(NA-k) ;
dP_k = @(z,k) nchoosek(NA,k) .* ps(z).^(k-1) .* (1-ps(z)).^(NA-k-1) .* dps(z) ;
Fbar_Ma_i{u,i} = @(z_) arrayfun(@(z) 0, z_) ;
% f_Ma_i{u,i} = Fbar_Ma_i{u,i} ;
for k = 0:NA
Fbar_Ma_i{u,i} = @(z_) arrayfun(@(z) Fbar_Ma_i{u,i}(z) + nchoosek(NA,k) .* (ps(z)).^k .* (1-ps(z)).^(NA-k) .* (Ups_i(z)).^k, z_) ;
F_Ma_i{u,i} = @(z_) arrayfun(@(z) 1-Fbar_Ma_i{u,i}(z), z_) ;
% f_Ma_i{u,i} = @(z_) arrayfun(@(z) f_Ma_i{u,i}(z) + dP_k(z,k) .* Ups_i(z) + P_k(z,k) .* dUps_i(z) , z_) ;
end
syms z
f_Ma_i{u,i} = eval( ['@(z)' char( diff( Ups_i(z) ) )] )
end
F_Ma{u,M} = @(z_) arrayfun(@(z) 1- ( min(realmax,Fbar_Ma_i{u,1}(z)) .* IND{u,1}(z) + min(realmax,Fbar_Ma_i{u,2}(z)) .* IND{u,2}(z) ...
+ min(realmax,Fbar_Ma_i{u,3}(z)) .* IND{u,3}(z) ), z_) ;
f_Ma{u,M} = @(z_) arrayfun(@(z) - ( min(realmax,f_Ma_i{u,1}(z)) .* IND{u,1}(z) + min(realmax,f_Ma_i{u,2}(z)) .* IND{u,2}(z) ...
+ min(realmax,f_Ma_i{u,3}(z)) .* IND{u,3}(z) ), z_) ;
end
end
The problem is that either when I computed manually and wrote it on MATLAB, when I used the symbolic variable, or when I tried to compute the derivative for each one of the pieces, I got errors like 'Index in position 1 exceeds array bounds (must not exceed 2)' or 'A and B must be floating-point scalar' for the function integral (although I'm using arrayfun).
Please let me know any suggested approach to compute this derivative or fix the errors I got.
Thank you in advance!

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by