Solving a nonlinear equation with numerical integration

Hello everybody,

I am pretty new in Matlab and struggling at the moment. I want to calculate two unknown parameters (x1 and x2) from two equations with integrals inside.

The Values P1, P2, S1, S2, C1, C2 are known and have the datatype double.

B(lambda) and e(lambda) are datasets (x,y) stored in arrays that describe the filter curves (also double).

Many thanks for your help.

Answer 1

[Torsten]

First method:

1. Fit functions B1(lambda), B2(lambda), e1(lambda) and e2(lambda) to your data in the different ranges (500 to 550 and 600 to 650).

2. Call fsolve as

fun = @(x)[P1-S1*integral(@(lambda)B1(lambda).*e1(lambda).*(1-exp(-x(1)./lambda.^1.39)).*c1./(lambda.^5.*(exp(c2./(lambda*x(2)))-1)),500,550),P2-S2*integral(@(lambda)B2(lambda).*e2(lambda).*(1-exp(-x(1)./lambda.^1.39)).*c1./(lambda.^5.*(exp(c2./(lambda*x(2)))-1)),600,650)];
x0 = [1, 1];
xsol = fsolve(fun,x0)

Second method:

function main
S1 = ...;
S2 = ...;
P1 = ...;
P2 = ...;
l = ...;
e = ...;
B = ...;
x0 = [1, 1];
xsol = fsolve(@(x)fun(x,S1,S2,P1,P2,l,e,B),x0);

function res = fun(x,S1,S2,P1,P2,l,e,B)
interpe = @(lambda)interp1(l,e,lambda);
interpB = @(lambda)interp1(l,B,lambda);
funI = @(lambda)interpe(lambda).*interpB(lambda).*(1-exp(-x(1)./lambda.^1.39)).*c1./(lambda.^5.*(exp(c2./(lambda*x(2)))-1));
res = [P1-S1*integral(funI,500,550),P2-S2*integral(funI,600,650)];

Best wishes

Torsten.

[Michael Bond]

Thank you very much for the fast answer. After trying to compile I get the following error message:

Undefined function 'fsolve' for input arguments of type 'function_handle'.

    %Error in Berechnung_T_KL (line 85)
    xsol = fsolve(fun,x0);

This link describes my error, after typing ver in the command window the optimization toolbox is listed: "Optimization Toolbox Version 7.5 (R2016b) " although I cant open it through Apps-->Optimization (another error appears in the command window: Undefined function or variable 'optimtool')

[Torsten]

You will need a rootfinder.

If you can't make fsolve run, you will have to program a simple newton method on your own or get code from the file exchange.

Best wishes

Torsten.

[Torsten]

Another possibility is to use "fminsearch" instead of "fsolve" for your problem:

function main
S1 = ...;
S2 = ...;
P1 = ...;
P2 = ...;
l = ...;
e = ...;
B = ...;
x0 = [1, 1];
xsol = fminsearch(@(x)fun(x,S1,S2,P1,P2,l,e,B),x0);

function res = fun(x,S1,S2,P1,P2,l,e,B)
interpe = @(lambda)interp1(l,e,lambda);
interpB = @(lambda)interp1(l,B,lambda);
funI = @(lambda)interpe(lambda).*interpB(lambda).*(1-exp(-x(1)./lambda.^1.39)).*c1./(lambda.^5.*(exp(c2./(lambda*x(2)))-1));
res = (P1-S1*integral(funI,500,550))^2+(P2-S2*integral(funI,600,650))^2;

[Michael Bond]

Hello Torsten, Update: I've managed to get a licence for the optimization toolbox, so now, I can use fsolve and I think its the right solver for my problem.

Here the important part of my code:

P1=100;
P2=150; 
S1=1;
S2=2;
fun = @(x)[P1-S1*integral(@(B1_x)B1(B1_x).*E(E_x).*(1-exp(-x(1)./B1_x.^1.39)).*c1./(B1_x.^5.*(exp(c2./(B1_x*x(2)))-1)),500,550),P2-S2*integral(@(B2_x)B2(B2_x).*E(E_x).*(1-exp(-x(1)./B2_x.^1.39)).*c1./(B2_x.^5.*(exp(c2./(B2_x*x(2)))-1)),600,650)];
x0 = [1, 1];
xsol = fsolve(fun,x0);

B1, B2 and E are 1751x2 double arrays, B1_x, B2_x and E_x are 1751x1 double arrays (x-Values), B1_y, B2_y, E_y have the same size.

After trying to compile I'll get another Errors:

Subscript indices must either be real positive integers or logicals.

Error in solve_equation>@(B1_x)B1(B1_x).*E(E_x).*(1-exp(-x(1)./B1_x.^1.39)).*c1./(B1_x.^5.*(exp(c2./(B1_x*x(2)))-1))

Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

Error in
solve_equation>@(x)[P1-S1*integral(@(B1_x)B1(B1_x).*E(E_x).*(1-exp(-x(1)./B1_x.^1.39)).*c1./(B1_x.^5.*(exp(c2./(B1_x*x(2)))-1)),400,750),P2-S2*integral(@(B2_x)B2(B2_x).*E(E_x).*(1-exp(-x(1)./B2_x.^1.39)).*c1./(B2_x.^5.*(exp(c2./(B2_x*x(2)))-1)),400,750)]

Error in fsolve (line 230)
            fuser = feval(funfcn{3},x,varargin{:});

Error in solve_equation (line 84)
xsol = fsolve(fun,x0);

Caused by:
    Failure in initial objective function evaluation. FSOLVE cannot continue.

[Torsten]

Use the second method.

l,e and B are assumed as one-dimensional arrays of the same size where l covers the complete range of integration for lambda and e(i) and B(i) are e and B, evaluated at l(i).

Best wishes

Torsten.

[Michael Bond]

Hello again:) I forgot an index in the attached equation in my first post, it should be B1 in the first equation and B2 in the second. I have fixed this issue in the code below. After compiling there are no errors, but where can I find the calculated results? They are not in the workspace or in the command window.

function main
S1 = 1.2;
S2 = 1.7;
P1 = 100;
P2 = 150;
l = B1_x;
e = E;
  x0 = [1, 1];
  xsol = fsolve(@(x)fun(x,S1,S2,P1,P2,l,e,B),x0);
  end 
  function res = fun(x,S1,S2,P1,P2,l,e,B1,B2)
  interpe = @(lambda)interp1(l,e,lambda);
  interpB1 = @(lambda)interp1(l,B1,lambda);
  interpB2 = @(lambda)interp1(l,B2,lambda);
  funI = @(lambda)interpe(lambda).*interpB1(lambda).*(1-exp(-x(1)./lambda.^1.39)).*c1./(lambda.^5.*(exp(c2./(lambda*x(2)))-1));
  funII = @(lambda)interpe(lambda).*interpB2(lambda).*(1-exp(-x(1)./lambda.^1.39)).*c1./(lambda.^5.*(exp(c2./(lambda*x(2)))-1));
  res = [P1-S1*integral(funI,400,750),P2-S2*integral(funII,400,750)]
  end

Many thanks! Michael

[Torsten]

Delete the semicolon after

xsol = fsolve(@(x)fun(x,S1,S2,P1,P2,l,e,B),x0);

But note that the lists of variables transferred to "fun" don't agree:

B becomes B1,B2.

Read again what l, e and B have to contain in my notation.

Best wishes

Torsten.

[Michael Bond]

There is still no error and no answer. I'll try to explain my variables clearlier:

E, B1, B2 are 1751x2 double arrays (x,y Data); with x =400:0.2:750. I have extracted the x and y values to new one dimensional arrays (1751x1 double), they are named as B1_x, B1_y, B2_x, B2_y, E_x, E_y. In principle it doesn't matter which x-Variable will be assigned to "e" because they have all the same values. Also the integration limits can be the same (400-750).

 function main
  S1 = 1.2;
  S2 = 1.7;
  P1 = 100;
  P2 = 150;
  l = B1_x;
  e = E_y;
  B1 = B1_y;
  B2 = B2_y;

x0 = [1, 1];
xsol = fsolve(@(x)fun(x,S1,S2,P1,P2,l,e,B1,B2),x0)
end 
function res = fun(x,S1,S2,P1,P2,l,e,B1,B2);
interpe = @(lambda)interp1(l,e,lambda);
interpB1 = @(lambda)interp1(l,B1,lambda);
interpB2 = @(lambda)interp1(l,B2,lambda);
funI = @(lambda)interpe(lambda).*interpB1(lambda).*(1-exp(-x(1)./lambda.^1.39)).*c1./(lambda.^5.*(exp(c2./(lambda*x(2)))-1));
funII = @(lambda)interpe(lambda).*interpB2(lambda).*(1-exp(-x(1)./lambda.^1.39)).*c1./(lambda.^5.*(exp(c2./(lambda*x(2)))-1));
res = [P1-S1*integral(funI,400,750),P2-S2*integral(funII,400,750)];
end

[Torsten]

Remove the two "end" statements in your code.

Remove the semicolon after

function res = fun(x,S1,S2,P1,P2,l,e,B1,B2);

Best wishes

Torsten.

[Michael Bond]

I've did it, it says "All functions in a script must be closed with an 'end'."