Working out a summation for fitting a curve

Hi, I want to work out a summation that goes as.

I want to fit this to a curve plotted as I on the y axis and t on the x axis using the curve fit app. My function for that goes as:

function I = CurveFitD(t, deltaQ, d)
%CurveFitD We want to fit the I(t)=f(t) to a custom function. 

%   For a bounded layer we have an expression which can yield us the
%   thickness of the polymer layer (in centimeters). We have to fit the transient to this
%   equation to calculate the layer thickness. 

D           = 1e-05; 

% deltaQ      = 4.2e-11;                                   % Couloumbs

I           = zeros(size(t));
for i    = 1:length(t)
    for k = 1:5
        I(i) = I(i) + (-(1)^(k)).*...
            (exp((-(k^2).*(d^2))))./(D*t(i));
    end
    I(i) = sqrt(1./t(i)).*(1+2.*I(i));
end

I        = sqrt(D/(pi))*(deltaQ/d).*I;

end

[Torsten]

I guess you have measurement data for I at certain times and you want to fit D ?

[Walter Roberson]

D_num = sym(1e-05)

D_num = 

delta_q_num = sym(4.2e-11)

delta_q_num = 

syms D d delta_q t
syms n integer
Pi = sym(pi);

inner = (-1).^n .* exp(-n.^2 .* d.^2 / (D .* t))

inner = 

inner_sum = symsum(inner, n, 1, inf)

inner_sum = 

multiplier = sqrt(D ./ (Pi * t)) .* delta_q / d

multiplier = 

I = multiplier * (1 + 2 * inner_sum)

I = 

Is = simplify(I)

Is = 

I_num = subs(Is, {D, delta_q}, {D_num, delta_q_num} )

I_num = 

We would need to know the value of d, and the range of time, in order to plot it.

It is not obvious in the above display, but there is a negative sign in the term.

[Hashim]

Hi, I have an array of I (the experimental output) vis a vis t (time) which I am trying to fit to this equation to find d. This function is being fed to the curve fitting app to find that d. For now I just want to know if the implementation of the equation in my function is correct or not. I also like that I am able to put in estimates for d and deltaQ. @Torsten has the right idea.

[Bruno Luong]

IMO your code is not correct

(exp((-(k^2).*(d^2))))./(D*t(i))

The denominator is outside exp argument, not as your math expression.

[Torsten]

Yes,

I(i) = I(i) + (-1).^k.*exp(-k.^2.*d.^2./(D*t(i)))   

instead of

I(i) = I(i) + (-(1)^(k)).*...
    (exp((-(k^2).*(d^2))))./(D*t(i));

[Torsten]

Bruno and I corrected your mistakes in the calculation of I(i).

Curious why you still make the same mistakes as in your first attempt to code the infinite series.

Worse:

Now you write

I(k) = I(k) + ...

instead of the correct

I(i) = I(i) + ...

[Hashim]

Corrected:

function I = CurveFitD(t, d)

D           = 1e-05; 

deltaQ      = 4.7e-6;                                   % Couloumbs

I           = zeros(size(t));
% Bounded Cottrell
for i    = 1:length(t)
    for k = 1:5
        I(i) = I(i) + (-1).^k.*exp(-k.^2.*d.^2./(D*t(i))); 
    end
    I(i) = sqrt(1./t(i)).*(1+2.*I(i));
end

I        = sqrt(D/(pi))*(deltaQ/d).*I;

end

[Bruno Luong]

It looks fine to me. You can also use MATLAB vectorization instead of for-loop:

D = 1e-05; 
deltaQ = 4.7e-6;
nmax = 5;
n = (1:nmax)';
t = reshape(t,1,[]);
I = sqrt(D./(pi*t)).*deltaQ/d.*(1 + 2*sum((-1).^n.*exp(-(n*d).^2./(D*t)),1));

[Hashim]

Thanks everyone!