How to fit with an infinite serie function ?

Hi all,

i am trying fit a complicated function involving an infinite sum. It looks like,

where M(t)/Meq is 'y', t is x. the variables are D and L, which are coefficient and thickness respectively. Is there anyway to do this in matlab?

I tried to write a script, but unfortunately it doesnot work for me.

%  FID fitting procedure
plot(Time,Mass,'bs','markersize',6); hold on;   %PLOT_1
xlabel('Time (m)','interpreter','latex','Fontsize',18); hold on;
ylabel ('Mt/Minf','interpreter','latex','Fontsize',18);hold on;

%Given parameters

coefficient = 0;         coeff_lowerlimit = 0;   coeff_upperlimit = 50;                            
thickness =  100;       thickness_lower = 0;      thickness_upper = 2000;

numIteration = 100;

d1 = coefficient;     
d2 = thickness;   

x1 = Time;       %assign x experimental data
y1 = Mass;    %assign y experimental data

%syms n
%f1 = 8/((2*n+1)^2*pi^2)*exp((-a1*x*(2*n+1)^2*pi^2./(4*a2^2));

%f2 = symsum(f1,n,1,inf);
syms n
data1 = 8/((2*n+1)^2*pi^2)*exp(-d1*x1*(2*n+1)^2*pi^2./(4*d2^2));
data2 = symsum (data1,n,1,inf);

%plot(TimeBase(START1:END1), data1,'r','LineWidth',3.0); hold on;                        %PLOT_2

d1l = coeff_lowerlimit+0.000001;
d1u = coeff_upperlimit+0.000001;
d2l = thickness_lower+0.000001;
d2u = thickness_upper+0.000001;

lb1 = [d1l d2l]; %lower limit for parameters matrix
ub1 = [d1u d2u]; %upper limit for parameters matrix
d01 = [d1 d2]; %initialise parameters matrix

options = optimset('Display', 'iter', 'TolFun',1e-25, 'TolX', 1e-25, 'MaxFunEvals', 1e10, 'MaxIter', numIteration);
[asym, resnorm, residual] = lsqcurvefit(@igsorp2,d01,x1,y1,lb1,ub1,options); %fitting procedure
display (asym);

d1sym = asym(1); d2sym = asym(2);   %best fit parameters matrix
syms n
data1sym = 8/((2*n+1)^2*pi^2)*exp(-d1sym*x1*(2*n+1)^2*pi^2./(4*d2sym^2));
data2sym = symsum (data1sym,n,1,inf);

plot(time, data2sym,'r','LineWidth',3.0); hold on;                      %PLOT_3

coefficient = d1sym;     %recalculate the parameters here
thickness = d2sym;
display (coefficient);
display (thickness);

Discrepancy1 = sqrt(sum(residual.*residual));
display (Discrepancy1);

Many thanks for your help!!!

Answer 1

[Roger Stafford]

If you substitute

   s = exp(-D*pi^2*t/(4*L^2))

your function becomes the sum of the infinite series

   f(s) = 1 - 8/pi^2*(s^1/1 + s^9/9 + s^25/25 + s^49/49 + ...)

and in turn it can be shown that this is equal to (<-- not true, see comment)

   1 - 8/pi^2*integral from 0 to s of 1/(2*u)*log((1+u)/(1-u)) with respect to u.

In other words, your infinite sum, which doesn't converge very rapidly, can be replaced by the above integral which should be much easier to compute.

(Correction: replaced 8/pi^3 with 8/pi^2)

(2nd Correction: The above integral is incorrect. See my comment below.)

[eric]

Do you mean use the substitution instead of using 'syms'? e.g.

      % code
d1sym = asym(1); d2sym = asym(2);   %best fit parameters matrix
ssym = exp(-d1sym*pi^2*x/(4*d2sym^2));
data1sym=1 - 8/pi^2*(ssym^1/1 + sssym^9/9 + s^25/25 + sssym^49/49 + ...)

    end

i guess i still not catch your point :(.

Thanks for your help.

[Roger Stafford]

The point of making the substitution

   s = exp(-D*pi^2*t/(4*L^2))

is that you need not deal with every possible combination of the variables, D, L, and t, but as far as computing 'f' is concerned you can deal with the single variable s.

The point of using the integral form for 'f' is that it should make evaluating it much easier than having to sum a rather slowly converging infinite series.

[eric]

Thanks for your help.

[Kuang-Sheng Chung]

Hi Eric, I face exact the same diffusion equation and want to estimate Msat, D, and L.

However, I used another approach. I assumed the diffusion process is simply an exponential when the absorption is over 50%.

I wonder that the method of Roger finally works or not. Do you have any further result?

[Roger Stafford]

My apologies to Eric, Antonio, and Kuang-Sheng. I have just discovered that my answer to the question posed here back in 2014 was wrong. As Antonio pointed out in a comment in 2015 which I have only just now read, the integral I suggested has a power series expansion of:

   s^1/1^2 + s^3/3^2 + s^5/5^2 + s^7/7^2 + ...

rather than the given series

   s^(1^2)/(1^2) + s^(3^2)/(3^2) + s^(5^2)/(5^2) + s^(7^2)/(7^2) + ...

and these are definitely not equivalent expressions. To my chagrin, somehow the square of the power disappeared during my manipulations at that time.

I do not know of a closed expression for this second, correct series, so presumably it would have to be evaluated by a direct, approximating finite summation. There is still an advantage for Eric in doing a fit varying only the single s value rather than D, t, and L, but the integral I suggested is incorrect. Again my apologies.

Answer 2

[Antonio Aguilera Miguel]

Hi Roger,

Thanks a lot for your answer to Eric. It is really useful for me but I am not able to find the same equality. I would be really grateful if you could correct me or tell me what is wrong!!

Turning back the integral from 0 to s of (1/(2*u))*log((1+u)/(1-u)) with respect u, and using an online integrator, the integral you proposed is equivalent to:

(1/2)*[-Li_2(-u)+Li_2(u)+ log(u)*(log((1-u)/(u+1))+log((u+1)/(1-u)))] being Li_2 Jonquière's functions.

Developing the down part, I got that log(u)*(log((1-u)/(u+1))+log((u+1)/(1-u))=0 so that means that the integral you proposed would result only (1/2)*[-Li_2(-u)+Li_2(u)].

Bearing in mind that Li_n(z)=sum(from k=0 to infinite)of (z^k)/(k^n); what it the definition of a Jonquière's function, I can develop these two terms in series and I got that:

[-Li_2(-u)+Li_2(u)]= 2*[u/1+(u^3)/(9)+(u^5)/(25)+(u^7)/(49)+....

So at the end, once I have applied the integral limits I got that the infinite sum would be in this way:

s/1+(s^3)/(9)+(s^5)/(25)+...

It still lacks me the numerator to elevate to the second power???

Could you or someone help me with this issue? I would be really grateful. Thanks you guys Wish you a very good day Antonio