Solving a nonlinear equation with numerical integration

2017 11 30
1 个回答
更新时间：2017 12 1
10 次查看（30 天）

0 个投票

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.

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2017-11-30
编辑：Torsten 2017-11-30

0 个投票

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.

10 个评论
显示 8更早的评论 隐藏 8更早的评论

Michael Bond
Michael Bond 2017-11-30
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
Torsten 2017-11-30
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
Torsten 2017-12-1
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
Michael Bond 2017-12-1
编辑：Michael Bond 2017-12-1
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
Torsten 2017-12-1
编辑：Torsten 2017-12-1
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
Michael Bond 2017-12-1
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
Torsten 2017-12-1
编辑：Torsten 2017-12-1
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
Michael Bond 2017-12-1
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
Torsten 2017-12-1
编辑：Torsten 2017-12-1
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
Michael Bond 2017-12-1
I've did it, it says "All functions in a script must be closed with an 'end'."

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Univariate Discrete Distributions 的更多信息

产品

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by