MATLAB: How to integrate a function with infinite series inside

infinite seriesnumerical integration

i have this function in the image that i want to integrate and it contains infinite series inside the intergral, i tried many method, but it did not work out.
any thought in how to do this intergration ?

Best Answer

  • 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; % 100 terms will be used in calculations

    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; % 100 terms will be used in calculations
    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)