data fitting by Integration with variable limit with a unknown parameter.

I want to fit experimental data with equation:

y = C(1) * 74.826 * (x/T(1))^3 * integration (t^4*exp(t)/(exp(t)-1)^2), 0, T(1)/x) + C(2) * 24.942 * (T(2)/x)^2 * exp(T(2)/x)/(exp(T(2)/x)-1)^2 + C(3) * 24.942 * (T(3)/x)^2 * *exp(T(3)/x)/(exp(T(3)/x)-1)^2

with six unknown parameters: C(1), T(1), C(2), T(2), C(3), T(3).

The integral part is giving error all times. How to write the programm please help.

Thank you.

Answer 1

[Bjorn Gustavsson]

You can do something like this:

y_fcn = @(CT,x) CT(1) * 74.826 * (x./CT(4)).^3 * integral(t.^4.*exp(t)./(exp(t)-1).^2), 0, CT(4)./x) + ...
                CT(2) * 24.942 * (CT(5)./x).^2  .* exp(CT(5)./x)./(exp(CT(5)./x)-1).^2 +  CT(3)*24.942*(CT(6)./x).^2.*exp(CT(6)./x)./(exp(CT(6)./x)-1).^2;
err_fcn = @(CT,x,y) sum((y-y_fcn(CT,x)).^2);
CT0 = [1 2 3 4 5 7]; % Some sensible initial guess for C and T, that we concatenate together as CT = [C,T]
CT_best = fminsearch(@(CT) err_fcn(CT,x_obs,y_obs),CT0)

You most likely have to use element-wise operations in the equation for y (I might have missed one or two multiplictions or power-operations). Then it should be fine to use ordinary least-square fitting of the C and T parameters (provided you have enough data-points, and they have the same standard deviation - but I expect you to know this already).

HTH

[Bjorn Gustavsson]

OK, I realise that I goofed in the definition of y_fcn - twice, the first fix you'll have to do is to make the first input argument to integral into a function, something like this:

y_fcn = @(CT,x) CT(1) * 74.826 * (x./CT(4)).^3 * integral(@(t) t.^4.*exp(t)./(exp(t)-1).^2, 0, CT(4)./x) + ...
                CT(2) * 24.942 * (CT(5)./x).^2  .* exp(CT(5)./x)./(exp(CT(5)./x)-1).^2 + ...
                CT(3)*24.942*(CT(6)./x).^2.*exp(CT(6)./x)./(exp(CT(6)./x)-1).^2;

note the "@(t) expression" - that is a definition of a anonymous function/dynamic function. The next problem is that your upper bound is now an array of values - that will also cause an error - exactly how to fix that I do not know since it is not clear what you want that expression to return - a single value of the integral from 0 to T(1)/x(end) or an array of the integrals from 0 to T(1)/x(1), from 0 to T(1)/x(2), ... from 0 to T(1)/x(end)?

First you have to clarify that.

[Shobha Gondh]

integral will return : an array of the integrals from 0 to T(1)/x(1), from 0 to T(1)/x(2), ... from 0 to T(1)/x(end).

error:

Unrecognized function or variable 'x_obs'.

Error in test12>@(CT)err_fcn(CT,x_obs,y_obs) (line 40)

CT_best = fminsearch(@(CT) err_fcn(CT,x_obs,y_obs),CT0);

Error in fminsearch (line 200)

fv(:,1) = funfcn(x,varargin{:});

Error in test12 (line 40)

CT_best = fminsearch(@(CT) err_fcn(CT,x_obs,y_obs),CT0);

and guess values are CT0 : [CT1 = 1, CT4 = 300, CT2 = 5, CT3 = 150, CT5 = 9, CT = 250]. Is this x_obs is x? means data at x-axis?

[Bjorn Gustavsson]

OK, we have 2 steps to fix.

1, the integral function only returns a value to an integral between 2 points, corresponding to the functions trapz and sum. Those two functions have cummulative variants, cumtrapz and cumsum, unfortunately integral has no such built-in-variant. However, I had a problem similar to yours and wrapped the calculations into a simple cumintegral-function. It is attached, and you should be able to replace "integral" with "cumintegral" in the definition of y_fcn:

y_fcn = @(CT,x) CT(1) * 74.826 * (x./CT(4)).^3 .* cumintegral(@(t) t.^4.*exp(t)./(exp(t)-1).^2, 0, CT(4)./x) + ...
                CT(2) * 24.942 * (CT(5)./x).^2  .* exp(CT(5)./x)./(exp(CT(5)./x)-1).^2 + ...
                CT(3)*24.942*(CT(6)./x).^2.*exp(CT(6)./x)./(exp(CT(6)./x)-1).^2;

Note that we have to change the matrix-multipy before cumintegral to a element-multiply.

2, In the call to fminsearch we use the dynamic anonymous-function capability to define a name-less function that is only a variable of 1 variable, CT, the other two inputs to err_fcn are taken at their values at the time of the call, in my code-snippet I named them x_obs and y_obs, those should be your x and y:

CT_best = fminsearch(@(CT) err_fcn(CT,x,y),CT0)

These modifications will give you something that returns a result. You might have to check that your y_fcn-function is a function that does what you expect it to.

[Shobha Gondh]

The equation is:

P = P0 + 4RC*(T/C)^5* integration (0, (C/T)) (x^5/(exp(x)-1)*(1-exp(-x))) dx

P0, R and C are constants.

With given values of P and T, finding these constants.

program in matlab for this :

Tv = linspace(200, 300, 10).'; % vector of T values

Pv = Tv.^2; % vector of P values

P_mdl = @(P0, R, C, TT) arrayfun(@(T) P0 + 4*R.*C.*(T/C).^5 * integral(@(x) x.^5./(exp(x)-1).*(1-exp(-x)), 0, (C/T)), TT);

sol = lsqcurvefit(@(param, T) P_mdl(param(1), param(2), param(3), T), rand(3,1), Tv, Pv);

P0_sol = sol(1);

R_sol = sol(2);

C_sol = sol(3);

This is same variable uper limit integration case but what are these Tv and Pv ? please help me to write program in this lsqcurvefit.

That cumintegral is also not giving the results. cumintegral is coming as undefined variable.

[Bjorn Gustavsson]

Ah, I forgot to attach the cumintegral function. Here it is...