For such series, if the sum is convergent, it is useful to approximate it with finite terms instead of using infinite series. You can use sum() and integral() functions to implement this equation. See the following code to see how it can be done
xf = 1;
xe = 2;
xd = 3;
sum1_terms = @(tda_,n) exp(-4.*n.^2.*pi^2.*tda_);
sum2_terms = @(tda_,n) exp(-4.*n.^2.*pi^2.*tda_).*sin(n.*pi.*xf./xe).*cos(n.*pi.*xd.*xf./xe)./(n.*pi.*xf./xe);
sum1 = @(tda_,N) sum(sum1_terms(tda_,1:N));
sum2 = @(tda_,N) sum(sum2_terms(tda_,1:N));
dPwd = @(tda_,N) (1+2*sum1(tda_,N)).*(1+2*sum2(tda_,N));
Pwd = @(tda, N) integral(@(tda_) dPwd(tda_, N), 0, tda, 'ArrayValued', 1);
N_value = 100;
tda_value = linspace(0,0.01,100);
Pwd_values = zeros(size(tda_value));
for i=1:numel(tda_value)
Pwd_values(i) = Pwd(tda_value(i), N_value);
end
plot(tda_value, Pwd_values)
Following code is a faster version of the above code (about 3x faster), but requires a bit deeper understanding of MATLAB functions
xf = 1;
xe = 2;
xd = 3;
sum1_terms = @(tda_,n) exp(-4.*n.^2.*pi^2.*tda_);
sum2_terms = @(tda_,n) exp(-4.*n.^2.*pi^2.*tda_).*sin(n.*pi.*xf./xe).*cos(n.*pi.*xd.*xf./xe)./(n.*pi.*xf./xe);
sum1 = @(tda_,N) sum(sum1_terms(tda_,(1:N).'));
sum2 = @(tda_,N) sum(sum2_terms(tda_,(1:N).'));
dPwd = @(tda_,N) (1+2*sum1(tda_,N)).*(1+2*sum2(tda_,N));
Pwd = @(tda, N) integral(@(tda_) dPwd(tda_, N), 0, tda);
N_value = 100;
tda_value = linspace(0,0.01,100);
Pwd_values = zeros(size(tda_value));
for i=1:numel(tda_value)
Pwd_values(i) = Pwd(tda_value(i), N_value);
end
plot(tda_value, Pwd_values)
Best Answer