Solve non linear equations system

Hi, I'm a newbie, I have never used Matlab but I have to use it to solve a complex problem as a part of a modelisation project. I need to solve 4 non linear equation systems, each one is a 2 equations system with two variables (x and y) and t is a constant. I would like to get the results for x and y depending on t. I have some boundaries conditions, 6.3*10^10>x>0, 6.3*10^10>y>0 and 300>t>0. One of my 4 system for example:

0=(10^(-41)*(10^9)*(e^((1-2*0.2)*160)-1))/(1-2*0.2)-((e^(-160)-1)/(0.2*(1-0.2)))*(0.91-0.35*(((12.6*10^(11)-t*y)^2-(12.6*10^11-t*x-t*y)^2)/(12.6*10^11-t*x-t*y)^2)-25/(12.6*10^11-t*y)-(x+y)*2*t*0.35*(((12.6*10^11-t*y)^2)/(12.6*10^11-t*x-t*y)^3))+(e^(-0.2*160)-1)/0.2*(-0.91*t+0.35*(2*(t^2)*x*(12.6*10^11-t*x-t*y)+(t^3)*(x^2))/((12.6*10^11-t*x-t*y)^2)-(25*t)/(12.6*10^11-t*y))

0=(2.7*10^(-23)*(e^((1-2*0.2)*90)-1))/(1-2*0.2)-((e^(-90)-1)/(0.2*(1-0.2)))*(0.91-0.79*(((12.6*10^(11)-t*x)^2-(12.6*10^11-t*y-t*x)^2)/(12.6*10^11-t*y-t*x)^2)-42/(12.6*10^11-t*x)-(y+x)*2*t*0.79*(((12.6*10^11-t*x)^2)/(12.6*10^11-t*y-t*x)^3))+(e^(-0.2*90)-1)/0.2*(-0.91*t+0.79*(2*(t^2)*y*(12.6*10^11-t*y-t*x)+(t^3)*(y^2))/((12.6*10^11-t*y-t*x)^2)-(42*t)/(12.6*10^11-t*x))

Could you please help me to use Matlab in order to solve my problem? Thank you very much

Answer 1

[jgg]

You need to do a couple of steps. First, make a function which returns the equation you want optimized:

 function [F] = fun_t(a,t)
 x = a(1);
 y = a(2);
 e = exp(1);
 F = zeros(2,1);
 F(1) = (10^(-41)*(10^9)*(e^((1-2*0.2)*160)-1))/(1-2*0.2)-((e^(-160)-1)/(0.2*(1-0.2)))*(0.91-0.35*(((12.6*10^(11)-t*y)^2-(12.6*10^11-t*x-t*y)^2)/(12.6*10^11-t*x-t*y)^2)-25/(12.6*10^11-t*y)-(x+y)*2*t*0.35*(((12.6*10^11-t*y)^2)/(12.6*10^11-t*x-t*y)^3))+(e^(-0.2*160)-1)/0.2*(-0.91*t+0.35*(2*(t^2)*x*(12.6*10^11-t*x-t*y)+(t^3)*(x^2))/((12.6*10^11-t*x-t*y)^2)-(25*t)/(12.6*10^11-t*y));
 F(2) = (2.7*10^(-23)*(e^((1-2*0.2)*90)-1))/(1-2*0.2)-((e^(-90)-1)/(0.2*(1-0.2)))*(0.91-0.79*(((12.6*10^(11)-t*x)^2-(12.6*10^11-t*y-t*x)^2)/(12.6*10^11-t*y-t*x)^2)-42/(12.6*10^11-t*x)-(y+x)*2*t*0.79*(((12.6*10^11-t*x)^2)/(12.6*10^11-t*y-t*x)^3))+(e^(-0.2*90)-1)/0.2*(-0.91*t+0.79*(2*(t^2)*y*(12.6*10^11-t*y-t*x)+(t^3)*(y^2))/((12.6*10^11-t*y-t*x)^2)-(42*t)/(12.6*10^11-t*x));
 end

Save this in your working directory. This is for your example system, and I'm assuming you've written it down properly and everything, since I just copied it from your post. You'll want to check this is correct.

I would first do the simplest thing, which is to just try and solve this directly without worrying about the constraints, then checking afterwards if they are met. You can do this using fsolve:

 t = 5; %some value you want
 func = @(a)fun_t(a,t)
 x0 = [10,10]; % a guess for the solution 
 x = fsolve(func,x0);

I think you might have a mistake in your system, since when I run this it evaluates to be perfectly flat. You'd better check that it's correct, but this just means changing F(1) and F(2) in the above function.

If you need the constraints, you can do the following instead of fsolve.

 lb = [0,0];
 ub = [6.3*10^10,6.3*10^10];
 rng default
 x0 = [10,10];
 [x res] = lsqnonlin(func,x0,lb,ub)

See also here for other techniques you can try. Some are easier than others.

[M N]

Thank you very much. I am actually using the trial version of Matlab. I tried to copy, past and save the file in the editor under the file name "myfunt.m". Then I copy and past your second command in the command window, but I still have a big red message when I try to run it. I think maybe something is wrong with the software.

Error: File: root2d.m Line: 10 Column: 1
This statement is not inside any function.
 (It follows the END that terminates the definition of the function "fun_t".)

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

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

Can't reload 'C:\Program Files\MATLAB\R2015b\bin\win64\mm_rmidd_mi.dll'

Error in rmidata.RmiData

Error in rmiml.RmiMlData

Error in rmiml.RmiMlData.getInstance

Error in rmiml.getAll

Can't reload 'C:\Program Files\MATLAB\R2015b\bin\win64\mm_rmidd_mi.dll'

Error in rmidata.RmiData

Error in rmiml.RmiMlData

Error in rmiml.RmiMlData.getInstance

Error in rmiml.visibleInToolstrip

Can't reload 'C:\Program Files\MATLAB\R2015b\bin\win64\mm_rmidd_mi.dll'

Error in rmidata.RmiData

Error in rmiml.RmiMlData

Error in rmiml.RmiMlData.getInstance

Error in rmiml.visibleInToolstrip

Can't reload 'C:\Program Files\MATLAB\R2015b\bin\win64\mm_rmidd_mi.dll'

Error in rmidata.RmiData

Error in rmiml.RmiMlData

Error in rmiml.RmiMlData.getInstance

Error in rmiml.visibleInToolstrip

Other try: I fixed the equation by replacing "e^" by "exp", because it's the exponential function.

Here is the script:

function [F] = fun_t(a,t)

   x = a(1);

   y = a(2);

   F = zeros(2,1);

   F(1) = (4.9*(10^9)*(exp((1-2*0.6)*100)-1))/(1-2*0.6)-((exp(-100)-1)/(0.6*(1-0.6)))*(0.91-0.35*(((12.6*10^(11)-t*y)^2-(12.6*10^11-t*x-t*y)^2)/(12.6*10^11-t*x-t*y)^2)-25/(12.6*10^11-t*y)-(x+y)*2*t*0.35*(((12.6*10^11-t*y)^2)/(12.6*10^11-t*x-t*y)^3))+(exp(-0.6*100)-1)/0.6*(-0.91*t+0.35*(2*(t^2)*x*(12.6*10^11-t*x-t*y)+(t^3)*(x^2))/((12.6*10^11-t*x-t*y)^2)-(25*t)/(12.6*10^11-t*y));

   F(2) = (330400*(exp((1-2*0.6)*50)-1))/(1-2*0.6)-((exp(-50)-1)/(0.6*(1-0.6)))*(0.91-0.79*(((12.6*10^(11)-t*x)^2-(12.6*10^11-t*y-t*x)^2)/(12.6*10^11-t*y-t*x)^2)-42/(12.6*10^11-t*x)-(y+x)*2*t*0.79*(((12.6*10^11-t*x)^2)/(12.6*10^11-t*y-t*x)^3))+(exp(-0.6*50)-1)/0.6*(-0.91*t+0.79*(2*(t^2)*y*(12.6*10^11-t*y-t*x)+(t^3)*(y^2))/((12.6*10^11-t*y-t*x)^2)-(42*t)/(12.6*10^11-t*x));

end

Trial>> t = 5; %some value you want
 func = @(a)fun_t(a,t)
 x0 = [10,10]; % a guess for the solution 
 x = fsolve(fun,x0);

func = 

      @(a)fun_t(a,t)

Not enough input arguments.

Error in fun (line 2)
 x = a(1);

[M N]

Would someone know how to fix it?

[jgg]

This works fine on my PC. Check to make sure you downloaded the file correctly. I've attached a version you can use.

[M N]

I tried again today, I also simplified the functions, so now it is :

function [F] = fun_t(a,t)
 x = a(1);
 y = a(2);
 F = zeros(2,1);
 F(1) = 8.20576*(12.6*10^11-t*y)*(12.6*10^11-t*x-t*y)^2+6.25*(((12.6*10^11-t*x-t*y)^2)*(0.91*(12.6*10^11-t*y)-25)-0.35*(12.6*10^11-t*y)*(2*12.6*10^11*y-2*(t^2)*x*y-(t^2)*(y^2)+t*(x+y)*(12.6*10^11-t*y)^2))-5*(0.35*(12.6*10^11-t*y)*((t^3)*(x^2)+2*(t^2)*x*(12.6*10^11-t*x-t*y))-t*((12.6*10^11-t*x-t*y)^2)*(25+0.91*(12.6*10^11-t*y)));
 F(2) = 12.7384*(12.6*10^11-t*x)*(12.6*10^11-t*x-t*y)^2+6.25*(((12.6*10^11-t*x-t*y)^2)*(0.91*(12.6*10^11-t*x)-42)-0.79*(12.6*10^11-t*x)*(2*12.6*10^11*x-2*(t^2)*x*y-(t^2)*(x^2)+t*(x+y)*(12.6*10^11-t*x)^2))-5*(0.79*(12.6*10^11-t*x)*((t^3)*(y^2)+2*(t^2)*y*(12.6*10^11-t*x-t*y))-t*((12.6*10^11-t*x-t*y)^2)*(42+0.91*(12.6*10^11-t*x)));
end

And then I get :

Trial>> t = 5; %some value you want
 func = @(a)fun_t(a,t)
 x0 = [10,10]; % a guess for the solution 
 x = fsolve(func,x0);

func = 

      @(a)fun_t(a,t)

No solution found.

fsolve stopped because the relative size of the current step is less than the
default value of the step size tolerance squared, but the vector of function values
is not near zero as measured by the default value of the function tolerance.

<stopping criteria details>

I don't really understand the meaning of the answer, is there anyone who know what it's about?

I also tried to fix t=20, but it still doesn't work :(

      function [F] = fun_t(a,t)
       x = a(1);
       y = a(2);
       F = zeros(2,1);
       F(1) = 8.20576*(12.6*10^11-20*y)*(12.6*10^11-20*x-20*y)^2+6.25*(((12.6*10^11-20*x-20*y)^2)*(0.91*(12.6*10^11-20*y)-25)-0.35*(12.6*10^11-20*y)*(2*12.6*10^11*y-2*(20^2)*x*y-(20^2)*(y^2)+20*(x+y)*(12.6*10^11-20*y)^2))-5*(0.35*(12.6*10^11-20*y)*((20^3)*(x^2)+2*(20^2)*x*(12.6*10^11-20*x-20*y))-20*((12.6*10^11-20*x-20*y)^2)*(25+0.91*(12.6*10^11-20*y)));
       F(2) = 12.7384*(12.6*10^11-20*x)*(12.6*10^11-20*x-20*y)^2+6.25*(((12.6*10^11-20*x-20*y)^2)*(0.91*(12.6*10^11-20*x)-42)-0.79*(12.6*10^11-20*x)*(2*12.6*10^11*x-2*(20^2)*x*y-(20^2)*(x^2)+20*(x+y)*(12.6*10^11-20*x)^2))-5*(0.79*(12.6*10^11-20*x)*((20^3)*(y^2)+2*(20^2)*y*(12.6*10^11-20*x-20*y))-20*((12.6*10^11-20*x-20*y)^2)*(42+0.91*(12.6*10^11-20*x)));
      end

      Trial>> t = 5; %some value you want
       func = @(a)fun_t(a,t)
       x0 = [10^7,10^7]; % a guess for the solution 
       x = fsolve(func,x0);

      func = 

          @(a)fun_t(a,t)

      Solver stopped prematurely.

      fsolve stopped because it exceeded the function evaluation limit,
      options.MaxFunEvals = 200 (the default value).

[jgg]

This means it can't find a solution. Nonlinear equations may not have solutions, or they can be difficult to find using the starting point. You need to ask yourself some questions:

• Do you know your system has a solution? It might not.
• Do you have an idea of where the solution might lie? You should try a starting point close to that if you can.
• You can also try several starting points, or one of the other algorithms supplied above.
• Are you sure your system is written down properly? You have terms like 8.20576*(12.6*10^11-t*y)*(12.6*10^11-t*x-t*y)^2 which is like 10^34 or more in magnitude which seems absurd.
• Can you simplify the system of equations at all (reducing the scale, for instance)?

Your Matlab code is working fine; you now just have to make sure the math problem you're trying to solve with it is going to work out.

[M N]

Well, x and y are extraction rate for Qatar and Iran over a common gas resource, their unit is MMBTU(gas quantity)/year, and the gas reserve is up to 12.6*10^11 MMBTU so it's ok to have such high values. I tried with xo=[10^7 10^7] but it's not working either.

[jgg]

This is a math problem you're going to have to solve; does this even have a solution? Your Matlab code is fine, you're just going to have to work on it for a while.

You should accept this answer so other people can see what you did to try and solve this problem.

Answer 2

[Krisda]

function [F] = fun_t(a,t)

x = a(1);

y = a(2);

e = exp(1);

F = zeros(2,1);

F(1) = (10^(-41)*(10^9)*(e^((1-2*0.2)*160)-1))/(1-2*0.2)-((e^(-160)-1)/(0.2*(1-0.2)))*(0.91-0.35*(((12.6*10^(11)-t*y)^2-(12.6*10^11-t*x-t*y)^2)/(12.6*10^11-t*x-t*y)^2)-25/(12.6*10^11-t*y)-(x+y)*2*t*0.35*(((12.6*10^11-t*y)^2)/(12.6*10^11-t*x-t*y)^3))+(e^(-0.2*160)-1)/0.2*(-0.91*t+0.35*(2*(t^2)*x*(12.6*10^11-t*x-t*y)+(t^3)*(x^2))/((12.6*10^11-t*x-t*y)^2)-(25*t)/(12.6*10^11-t*y));

F(2) = (2.7*10^(-23)*(e^((1-2*0.2)*90)-1))/(1-2*0.2)-((e^(-90)-1)/(0.2*(1-0.2)))*(0.91-0.79*(((12.6*10^(11)-t*x)^2-(12.6*10^11-t*y-t*x)^2)/(12.6*10^11-t*y-t*x)^2)-42/(12.6*10^11-t*x)-(y+x)*2*t*0.79*(((12.6*10^11-t*x)^2)/(12.6*10^11-t*y-t*x)^3))+(e^(-0.2*90)-1)/0.2*(-0.91*t+0.79*(2*(t^2)*y*(12.6*10^11-t*y-t*x)+(t^3)*(y^2))/((12.6*10^11-t*y-t*x)^2)-(42*t)/(12.6*10^11-t*x));

end

t = 5; %some value you want

func = @(a)fun_t(a,t)

x0 = [10,10]; % a guess for the solution

x = fsolve(func,x0);

lb = [0,0];

ub = [6.3*10^10,6.3*10^10];

rng default

x0 = [10,10];

[x res] = lsqnonlin(func,x0,lb,ub)