### 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.todaysummary.com.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.todaysummary.com.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.todaysummary.com.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.todaysummary.com/q_matlab_59019.html*

All Comments

Leave a comment...

- 7 Comments
- In article <ef1468b.-1.matlab.todaysummary.com.webx.raydaftYaTP>,
Ilya <theweremonkey.matlab.todaysummary.com.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.todaysummary.com.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.todaysummary.com.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.todaysummary.com.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.todaysummary.com.(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.todaysummary.com.(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.todaysummary.com.(x) quad(.matlab.todaysummary.com.(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.todaysummary.com.(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.todaysummary.com.(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.todaysummary.com.(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

- In article <ef1468b.-1.matlab.todaysummary.com.webx.raydaftYaTP>,
- 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.todaysummary.com.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.todaysummary.com.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.todaysummary.com.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

- I am a novice MatLab user, so I have little or no insight re: a
- 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.todaysummary.com.(sol) quad(.matlab.todaysummary.com.(c,sol) 1./sqrt(c.^2 -

sol(1).^2),sol(1)*(1+eps),x.*(1-eps),1e-10,[],x(1));

G5 = .matlab.todaysummary.com.(sol) quad(.matlab.todaysummary.com.(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.todaysummary.com.(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

- Hello again,
- 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

- Oh, and I keep forgetting to mention the following constraints:
- 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.todaysummary.com.(sol) quad(.matlab.todaysummary.com.(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.todaysummary.com.(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.todaysummary.com.(sol) quad(.matlab.todaysummary.com.(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.todaysummary.com.(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.todaysummary.com.(sol) quad(.matlab.todaysummary.com.(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.todaysummary.com.(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.todaysummary.com.(sol) quad(.matlab.todaysummary.com.(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.todaysummary.com.(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.todaysummary.com.(sol) quad(.matlab.todaysummary.com.(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.todaysummary.com.(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.todaysummary.com.(sol) quad(.matlab.todaysummary.com.(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.todaysummary.com.(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

- sigh... in the quad calls I forgot to substitute x with sol, so I
- 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.todaysummary.com.(sol) quad(.matlab.todaysummary.com.(c,sol) 1./sqrt(c.^2 -

sol(1).^2),sol(1)*(1+eps),x(j).*(1-eps),1e-10,[],sol(1));

G5 = .matlab.todaysummary.com.(sol) quad(.matlab.todaysummary.com.(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.todaysummary.com.(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

- I had someone who wasnt passing out from exhaustion look at the
- 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

- Dear Ilya,