Tags: adefinite, contains, equations, fsolve, integral, matlab, nonlinear, programming, series, solve, system, wasplanning
system of nonlinear equations
On Programmer » Matlab
18,589 words with 7 Comments; publish: Tue, 06 May 2008 19:06:00 GMT; (20062.50, « »)
I need to solve a series of 2 nonlinear equations, which I was
planning on doing with fsolve. One of the equations contains a
definite integral, the integrand and lower limit of which contain one
of the 2 variables that I am trying to solve for. So I can't define a
function to use with quad, because that definition will contain an
undefined variable that I am using in fsolve just like it is shown in
the tutorial.
Just as an example of what I am looking for:
Lets say that
equ1 : x(1) + x(2)*G1 = 0
equ2 : (2/pi)*x(2)*(3+x(1))-1 = 0
where
G1 = quad(.matlab.itags.org.integ,x(1),0.5)
and
integ(c) = 1/sqrt(c^2 - x(1)^2)
How could I find the solution to this system of equations?
Right now I am trying to do something like this:
function F=funcs(x)
F = [x(1) + x(2) * quad(.matlab.itags.org.integ,x(1),0.5);
(2/pi)*x(2)*(3+x(1))-1]
function G1 = integ(c)
G1 = 1/sqrt(c^2 - x(1)^2);
end
end
Then at the command I type:
x0 = [2.571 2.571];
x = fsolve(.matlab.itags.org.funcs, x0);
I know that this doesn't work and I basically understand why. I just
need to know if there is a way to solve this system of nonlinear
equations and what this way is.
http://matlab.itags.org/q_matlab_59019.html
All Comments
Leave a comment...
- 7 Comments
- In article <ef1468b.-1.matlab.itags.org.webx.raydaftYaTP>,
Ilya <theweremonkey.matlab.itags.org.hotmail.com> wrote:
> I need to solve a series of 2 nonlinear equations, which I was
> planning on doing with fsolve. One of the equations contains a
> definite integral, the integrand and lower limit of which contain one
> of the 2 variables that I am trying to solve for. So I can't define a
> function to use with quad, because that definition will contain an
> undefined variable that I am using in fsolve just like it is shown in
> the tutorial.
> Just as an example of what I am looking for:
> Lets say that
> equ1 : x(1) + x(2)*G1 = 0
> equ2 : (2/pi)*x(2)*(3+x(1))-1 = 0
> where
> G1 = quad(.matlab.itags.org.integ,x(1),0.5)
> and
> integ(c) = 1/sqrt(c^2 - x(1)^2)
> How could I find the solution to this system of equations?
> Right now I am trying to do something like this:
> function F=funcs(x)
> F = [x(1) + x(2) * quad(.matlab.itags.org.integ,x(1),0.5);
> (2/pi)*x(2)*(3+x(1))-1]
> function G1 = integ(c)
> G1 = 1/sqrt(c^2 - x(1)^2);
> end
> end
> Then at the command I type:
> x0 = [2.571 2.571];
> x = fsolve(.matlab.itags.org.funcs, x0);
> I know that this doesn't work and I basically understand why. I just
> need to know if there is a way to solve this system of nonlinear
> equations and what this way is.
We should be able to do it all anonymously.
G1 is a function of x(1) and c.
K = .matlab.itags.org.(c,x_1) 1./sqrt(c.^2 - x_1.^2);
We can integrate it for fixed x_1 by
x_1 = .2;
quad(K,x_1,.5,[],[],x_1)
Here is the first problem. We get a divide by zero,
because of the singularity at c = x(1).
Warning: Divide by zero.
> In .matlab.itags.org.(c,x_1) 1./sqrt(c.^2 - x_1.^2)
In quad at 62
It still works however. Quad is able to survive the
singularity.
ans =
1.56680960164147
We can try to avoid it completely by offsetting the
lower limit by a tiny relative amount.
quad(K,x_1*(1+eps),.5,1e-10,[],x_1)
ans =
1.56679921547824
One problem to consider is the integral is only defined
for x(1) in the open interval (0,.5), so we will need to
constrain x_1. There are many ways of doing it. Simplest
in an fsolve context is to use lsqnonlin, which allows
bound constraints.
First, lets write G1 as a function of x. I'll make some
of this simple by putting K inline with respect to G1.
G1 = .matlab.itags.org.(x) quad(.matlab.itags.org.(c,x) 1./sqrt(c.^2 - x(1).^2),
x(1)*(1+eps),.5,1e-10,[],x(1))
Now lets do the optimization, using the bound constraints
to avoid singularities.
Xstart = [.25 1];
LB = [1000*eps,-inf];
UB = [.5 - 1000*eps,inf];
options = optimset('disp','iter');
OBJ = .matlab.itags.org.(x) [x(1) + x(2)*G1(x), (2/pi)*x(2)*(3+x(1))-1];
Xfinal = lsqnonlin(OBJ,Xstart,LB,UB,options)
This works for a while, until we get a NaN as the result
of the finite differencing. The tiny change for the finite
difference pushed a point past the constraint limits. Then
it gets upset.
We can try to fix this by the artful use of diffmaxchange,
and by making the offsets on the limits a bit wider.
Xstart = [.25 1];
LB = [1e-4,-inf];
UB = [.5 - 1e-4,inf];
options = optimset('disp','iter','diffmaxchange',1
.e-8);
OBJ = .matlab.itags.org.(x) [x(1) + x(2)*G1(x), (2/pi)*x(2)*(3+x(1))-1];
Xfinal = lsqnonlin(OBJ,Xstart,LB,UB,options)
This gives us a nice result. (Note that I always use
'display' set to 'iter' as long as I am testing my code.
I only change it to off or final when I am confident of
its behavior.)
Norm of First-order
Iteration Func-count f(x) step optimality
CG-iterations
0 3 3.59815 4.28
1 6 0.514167 0.686304 0.0489
1
2 9 0.393858 0.37421 0.465
1
3 12 0.277751 0.14719 0.109
1
4 15 0.260506 0.0331703 0.0212
1
5 18 0.25897 0.00721908 0.00206
1
6 21 0.258935 0.000965385 4.29e-05
1
7 24 0.258935 4.05654e-05 8.09e-08
1
Optimization terminated: first-order optimality less than OPTIONS.TolFun,
and no negative/zero curvature detected in trust region model.
Xfinal =
0.499899999999933 0.4467616718016
We could have used other optimizers. One simple
solution that can avoid the diffmaxchange problem
is to use fminsearchbnd, from matlab central. We
still need to worry a little bit, because
fminsearchbnd is written to allow the boundary
values as valid parameter values.
Xstart = [.25 1];
LB = [1e-8,-inf];
UB = [.5 - 1e-8,inf];
options = optimset('disp','iter','tolfun',1.e-10);
OBJ = .matlab.itags.org.(x) sum([x(1) + x(2)*G1(x), (2/pi)*x(2)*(3+x(1))-1].^2);
Xfinal = fminsearchbnd(OBJ,Xstart,LB,UB,options)
After a few more iterations than lsqnonlin required,
fminsearchbnd also succeeds.
Iteration Func-count min f(x) Procedure
0 1 3.59815
1 3 3.59752 initial simplex
2 5 2.80246 expand
3 7 2.44918 expand
(... snip ...)
88 169 0.25009 contract inside
89 171 0.25009 contract outside
90 173 0.25009 contract inside
91 175 0.25009 contract inside
Optimization terminated:
the current x satisfies the termination criteria using OPTIONS.TolX of
1.000000e-04
and F(X) satisfies the convergence criteria using OPTIONS.TolFun of
1.000000e-10
Xfinal =
0.499999989999995 0.448779723638837
Since I allowed fminseachbnd to approach its limits
more tightly than I did lsqnonlin, we actually got a
slightly better result.
HTH,
John D'Errico
The best material model of a cat is another, or
preferably the same, cat.
A. Rosenblueth, Philosophy of Science, 1945
#1; Tue, 06 May 2008 19:07:00 GMT
- I am a novice MatLab user, so I have little or no insight re: a
tolerable numerical solution. However, there is something you might
want to try - to lose the integral term G1 in your system (I'm
currently by the lower limit of integration), parameterize
x(1) like:
x(1) = c*sin(u) so that
... dx(1) = c*cos(u)*du
G1 then turns into arcsin(x(1)/c),
and then evaluate at the appropriate limits.
However if there are poles to integrate through, you'll have to turn
G1 into a Cauchy integral. Regardless, this is a complicated system
to solve. Any nonlinear term + Murphy's Law will guarantee such.
Good luck!
Ilya wrote:
>
> I need to solve a series of 2 nonlinear equations, which I was
> planning on doing with fsolve. One of the equations contains a
> definite integral, the integrand and lower limit of which contain
> one
> of the 2 variables that I am trying to solve for. So I can't define
> a
> function to use with quad, because that definition will contain an
> undefined variable that I am using in fsolve just like it is shown
> in
> the tutorial.
> Just as an example of what I am looking for:
> Lets say that
> equ1 : x(1) + x(2)*G1 = 0
> equ2 : (2/pi)*x(2)*(3+x(1))-1 = 0
> where
> G1 = quad(.matlab.itags.org.integ,x(1),0.5)
> and
> integ(c) = 1/sqrt(c^2 - x(1)^2)
> How could I find the solution to this system of equations?
> Right now I am trying to do something like this:
> function F=funcs(x)
> F = [x(1) + x(2) * quad(.matlab.itags.org.integ,x(1),0.5);
> (2/pi)*x(2)*(3+x(1))-1]
> function G1 = integ(c)
> G1 = 1/sqrt(c^2 - x(1)^2);
> end
> end
> Then at the command I type:
> x0 = [2.571 2.571];
> x = fsolve(.matlab.itags.org.funcs, x0);
> I know that this doesn't work and I basically understand why. I
> just
> need to know if there is a way to solve this system of nonlinear
> equations and what this way is.
#2; Tue, 06 May 2008 19:08:00 GMT
- Hello again,
Thank you very much John for all of your help. I tried using your
method, but unfortunately, I am not familiar with how it solves this
system of equations, so I am having trouble tuning it to what I need.
Below I have included the code that I made by putting the information
I need into your code. Below the code is the result that I am
getting. I know for a fact that I shouldn't be getting complex
numbers here and I have checked my equations quite a few times. What
exactly is Xstart and what are UB and LB the boundaries of?
---
CODE
---
clear;
omega = 0.7;
x = 0.6;
G3 = .matlab.itags.org.(sol) quad(.matlab.itags.org.(c,sol) 1./sqrt(c.^2 -
sol(1).^2),sol(1)*(1+eps),x.*(1-eps),1e-10,[],x(1));
G5 = .matlab.itags.org.(sol) quad(.matlab.itags.org.(c,sol) sqrt(c.^2 -
sol(1).^2),sol(1)*(1+eps),x.*(1-eps),1e-10,[],x(1));
Xstart = [0.25 1];
LB = [1e-4,-inf];
UB = [x - 1e-4,inf];
options = optimset('disp','iter','diffmaxchange',1
.e-8);
OBJ = .matlab.itags.org.(sol)
[(omega.*sol(1))-(sol(2).*x.^3).*(sol(2).*((2/pi).*asin(sol(1)./x)+(1/
pi).*(-(sol(1)./x) + ((2.*sol(1)./x).*G3(sol)) - (1/(sol(1).*x)) +
((2./(sol(1).*x)).*G5(sol)))) - 1).^3,
(2/pi).*sol(2).*(asin(sol(1)./x)-(1./sol(1)).*(x-sqrt(x.^2 -
sol(1)))) - 1];
Xfinal = lsqnonlin(OBJ,Xstart,LB,UB,options)
---
RESULTS
---
Norm of First-order
Iteration Func-count f(x) step optimality
CG-iterations
0 3 37.5861 127
1 6 5.1747 0.350163 14
1
2 9 1.62863 0.298574 1.92
1
3 12 1.01456 0.38452 0.463
1
4 15 0.595433 2.1263 1.18
1
5 18 0.595433 2.17781 1.18
1
6 21 0.595433 0.544453 1.18
0
7 24 0.208778 0.136115 0.528
0
8 27 0.165713 0.272227 0.301
1
9 30 0.099203 0.544453 0.255
1
10 33 0.099203 1.26643 0.255
1
11 36 0.0628192 0.272227 0.0661
0
12 39 0.0183669 0.544453 0.367
1
13 42 0.000801249 0.458164 0.223
1
14 45 8.84871e-007 0.00999008 0.00731
1
15 48 1.40196e-012 0.00123489 9.73e-006
1
Optimization terminated: relative function value
changing by less than OPTIONS.TolFun.
Xfinal =
0.3364 - 0.2822i -3.3073 + 1.7069i
---
God I wish these windows were wider than they are.
#3; Tue, 06 May 2008 19:09:00 GMT
- Oh, and I keep forgetting to mention the following constraints:
0 <= sol(1) <= x
and
sol(2) >= pi/(pi-2)
#4; Tue, 06 May 2008 19:10:00 GMT
- sigh... in the quad calls I forgot to substitute x with sol, so I
have new incorrect results. Also, omega will always be equal to 1,
not 0.7 as I have in the code above.
Norm of First-order
Iteration Func-count f(x) step optimality
CG-iterations
0 3 6.85862 12.1
1 6 1.8807 0.506582 1.68
1
2 9 0.94928 0.592411 0.434
1
3 12 0.94928 2.18119 0.434
1
4 15 0.538244 0.545297 0.264
0
5 18 0.201212 1.09059 0.137
1
6 21 0.201212 1.33982 0.137
1
7 24 0.134853 0.334954 0.116
0
8 27 0.134853 0.775057 0.116
1
Warning: Minimum step size reached; singularity possible.
> In quad at 88
In test>.matlab.itags.org.(sol) quad(.matlab.itags.org.(c,sol) 1./sqrt(c.^2 -
sol(1).^2),sol(1)*(1+eps),x.*(1-eps),1e-10,[],sol(1)) at 5
In test>.matlab.itags.org.(sol)
[(omega.*sol(1))-(sol(2).*x.^3).*(sol(2).*((2/pi).*asin(sol(1)./x)+(1/
pi).*(-(sol(1)./x) + ((2.*sol(1)./x).*G3(sol)) - (1/(sol(1).*x)) +
((2./(sol(1).*x)).*G5(sol)))) - 1).^3,
(2/pi).*sol(2).*(asin(sol(1)./x)-(1./sol(1)).*(x-sqrt(x.^2 -
sol(1)))) - 1] at 12
In optim\private\snls at 395
In optim\private\lsqncommon at 222
In lsqnonlin at 147
In test at 14
Warning: Minimum step size reached; singularity possible.
> In quad at 88
In test>.matlab.itags.org.(sol) quad(.matlab.itags.org.(c,sol) 1./sqrt(c.^2 -
sol(1).^2),sol(1)*(1+eps),x.*(1-eps),1e-10,[],sol(1)) at 5
In test>.matlab.itags.org.(sol)
[(omega.*sol(1))-(sol(2).*x.^3).*(sol(2).*((2/pi).*asin(sol(1)./x)+(1/
pi).*(-(sol(1)./x) + ((2.*sol(1)./x).*G3(sol)) - (1/(sol(1).*x)) +
((2./(sol(1).*x)).*G5(sol)))) - 1).^3,
(2/pi).*sol(2).*(asin(sol(1)./x)-(1./sol(1)).*(x-sqrt(x.^2 -
sol(1)))) - 1] at 12
In sfdnls at 90
In optim\private\snls at 404
In optim\private\lsqncommon at 222
In lsqnonlin at 147
In test at 14
9 30 0.134853 0.167477 0.116
0
Warning: Minimum step size reached; singularity possible.
> In quad at 88
In test>.matlab.itags.org.(sol) quad(.matlab.itags.org.(c,sol) 1./sqrt(c.^2 -
sol(1).^2),sol(1)*(1+eps),x.*(1-eps),1e-10,[],sol(1)) at 5
In test>.matlab.itags.org.(sol)
[(omega.*sol(1))-(sol(2).*x.^3).*(sol(2).*((2/pi).*asin(sol(1)./x)+(1/
pi).*(-(sol(1)./x) + ((2.*sol(1)./x).*G3(sol)) - (1/(sol(1).*x)) +
((2./(sol(1).*x)).*G5(sol)))) - 1).^3,
(2/pi).*sol(2).*(asin(sol(1)./x)-(1./sol(1)).*(x-sqrt(x.^2 -
sol(1)))) - 1] at 12
In optim\private\snls at 395
In optim\private\lsqncommon at 222
In lsqnonlin at 147
In test at 14
Warning: Minimum step size reached; singularity possible.
> In quad at 88
In test>.matlab.itags.org.(sol) quad(.matlab.itags.org.(c,sol) 1./sqrt(c.^2 -
sol(1).^2),sol(1)*(1+eps),x.*(1-eps),1e-10,[],sol(1)) at 5
In test>.matlab.itags.org.(sol)
[(omega.*sol(1))-(sol(2).*x.^3).*(sol(2).*((2/pi).*asin(sol(1)./x)+(1/
pi).*(-(sol(1)./x) + ((2.*sol(1)./x).*G3(sol)) - (1/(sol(1).*x)) +
((2./(sol(1).*x)).*G5(sol)))) - 1).^3,
(2/pi).*sol(2).*(asin(sol(1)./x)-(1./sol(1)).*(x-sqrt(x.^2 -
sol(1)))) - 1] at 12
In sfdnls at 90
In optim\private\snls at 404
In optim\private\lsqncommon at 222
In lsqnonlin at 147
In test at 14
10 33 0.130226 0.0418692 0.407
0
Warning: Minimum step size reached; singularity possible.
> In quad at 88
In test>.matlab.itags.org.(sol) quad(.matlab.itags.org.(c,sol) 1./sqrt(c.^2 -
sol(1).^2),sol(1)*(1+eps),x.*(1-eps),1e-10,[],sol(1)) at 5
In test>.matlab.itags.org.(sol)
[(omega.*sol(1))-(sol(2).*x.^3).*(sol(2).*((2/pi).*asin(sol(1)./x)+(1/
pi).*(-(sol(1)./x) + ((2.*sol(1)./x).*G3(sol)) - (1/(sol(1).*x)) +
((2./(sol(1).*x)).*G5(sol)))) - 1).^3,
(2/pi).*sol(2).*(asin(sol(1)./x)-(1./sol(1)).*(x-sqrt(x.^2 -
sol(1)))) - 1] at 12
In optim\private\snls at 395
In optim\private\lsqncommon at 222
In lsqnonlin at 147
In test at 14
Warning: Minimum step size reached; singularity possible.
> In quad at 88
In test>.matlab.itags.org.(sol) quad(.matlab.itags.org.(c,sol) 1./sqrt(c.^2 -
sol(1).^2),sol(1)*(1+eps),x.*(1-eps),1e-10,[],sol(1)) at 5
In test>.matlab.itags.org.(sol)
[(omega.*sol(1))-(sol(2).*x.^3).*(sol(2).*((2/pi).*asin(sol(1)./x)+(1/
pi).*(-(sol(1)./x) + ((2.*sol(1)./x).*G3(sol)) - (1/(sol(1).*x)) +
((2./(sol(1).*x)).*G5(sol)))) - 1).^3,
(2/pi).*sol(2).*(asin(sol(1)./x)-(1./sol(1)).*(x-sqrt(x.^2 -
sol(1)))) - 1] at 12
In sfdnls at 90
In optim\private\snls at 404
In optim\private\lsqncommon at 222
In lsqnonlin at 147
In test at 14
11 36 0.127118 0.0418692 0.0556
1
12 39 0.127048 0.0837385 0.0329
1
13 42 0.126813 0.0209346 0.00675
1
14 45 0.126813 0.0418692 0.00675
1
15 48 0.126782 0.0104673 0.00282
0
16 51 0.126782 0.0104673 0.00282
1
17 54 0.126779 0.00261683 0.000975
0
18 57 0.126779 0.00523366 0.000975
1
19 60 0.126779 0.00130841 0.000745
0
Optimization terminated: relative function value
changing by less than OPTIONS.TolFun.
Xfinal =
0.3522 -1.8616
#5; Tue, 06 May 2008 19:11:00 GMT
- I had someone who wasnt passing out from exhaustion look at the
equations that I typed in and he found 2 errors. Now, I get almost
the right solutions. The values for sol(1) fall within the correct
range, but the higher the value of x gets, the farther away from zero
the solutions to the equations get. For x = 0.1 the value of equation
2 for the given equations is -4.8948e-005, but for x = 0.8, the
result is -0.0218. Now, as you can imaging, this is a bit too far
off. I am assuming that I am getting these results because I don't
know how to tune the parameters of lsqnonlin to my specific
situation.
Here is my code. The equations are correct, but I think I need to
alter either UB and LB, or Xstart and I dont know exactly what any of
them do. Help?
---
clear;
omega = 1;
x = 0.1:0.1:1;
for j=1:10
G3 = .matlab.itags.org.(sol) quad(.matlab.itags.org.(c,sol) 1./sqrt(c.^2 -
sol(1).^2),sol(1)*(1+eps),x(j).*(1-eps),1e-10,[],sol(1));
G5 = .matlab.itags.org.(sol) quad(.matlab.itags.org.(c,sol) sqrt(c.^2 -
sol(1).^2),sol(1)*(1+eps),x(j).*(1-eps),1e-10,[],sol(1));
Xstart = [0.25 1];
LB = [1e-4,-inf];
UB = [x(j) - 1e-4,inf];
options = optimset('disp','iter','diffmaxchange',1
.e-8);
OBJ = .matlab.itags.org.(sol)
[(omega.*sol(1))-(sol(2).*x(j).^3).*(sol(2).*((2/pi).*asin(sol(1)./x(j
))+(1/pi).*(-(sol(1)./x(j)) + ((2.*sol(1)./x(j)).*G3(sol)) -
((1./(sol(1).*x(j))).*(x(j).^2-sol(1).^2)) +
((2./(sol(1).*x(j))).*G5(sol)))) - 1).^3,
(2/pi).*sol(2).*(asin(sol(1)./x(j))-(1./sol(1)).*(x(j)-sqrt(x(j).^2 -
sol(1).^2))) - 1];
Xfinal = lsqnonlin(OBJ,Xstart,LB,UB,options)
R(j,:)=Xfinal;
end
---
#6; Tue, 06 May 2008 19:12:00 GMT
- Dear Ilya,
I have a problem that is similar (partitially) to yours. Could you
please give me some help with it? It would be great if you send me
your m-files.
Best regards
#7; Tue, 06 May 2008 19:13:00 GMT