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; (20078.13, « »)

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.

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
    • 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
    • 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
    • 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.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
    • 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
    • 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