help fsolve not getting correct answers

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

help fsolve not getting correct answers

Imane E
Hello,

I am pretty new to Octave and I need some help figuring out a code.
I wrote this code:
function y = f(x)
 
  y = zeros(2,1);
 
  P1=100000;
  nu1=.8333;
  gamma=1.4;
  mu2=(gamma-1)/(gamma+1);
  q=1200e3;
  E=2*q*mu2/(P1*nu1);
 
 
  a = 1 + x(1)^2/(P1*nu1)*(1-x(2));
  b =  -mu2 + (1 - mu2^2 + E)/(x(2) - mu2);
 
  c = -x(1)^2/(P1*nu1);
  d = -(1 - mu2^2 + E)/((x(2) - mu2)^2);
 
  y(1) = a - b ;
  y(2) = c - d ;
 
  endfunction

And then used:
[x, fval, info] = fsolve (@f, [1;2])
to solve for x(1) and x(2) when a = b and c = d (where c and d are respective derivatives of a and b).

I know that the x values are supposed to be x(1)=1000 and x(2)=0.5, as I previously solved this system by hand. However, fsolve gives me:
x =

  -73.3172
    9.6259

Is there anything I can do to improve he code so that it gives the correct answers, or something close?



Reply | Threaded
Open this post in threaded view
|

RE: help fsolve not getting correct answers

Tony Richardson
> Subject: help fsolve not getting correct answers
>
> Hello,
>
> I am pretty new to Octave and I need some help figuring out a code.
> I wrote this code:
> function y = f(x)

>  y = zeros(2,1);

>  P1=100000;
>  nu1=.8333;
>  gamma=1.4;
>  mu2=(gamma-1)/(gamma+1);
>  q=1200e3;
>  E=2*q*mu2/(P1*nu1);


>  a = 1 + x(1)^2/(P1*nu1)*(1-x(2));
>  b =  -mu2 + (1 - mu2^2 + E)/(x(2) - mu2);

>  c = -x(1)^2/(P1*nu1);
>  d = -(1 - mu2^2 + E)/((x(2) - mu2)^2);

>  y(1) = a - b ;
>  y(2) = c - d ;

>  endfunction
>
> And then used:
> [x, fval, info] = fsolve (@f, [1;2])
> to solve for x(1) and x(2) when a = b and c = d (where c and d are respective derivatives of a and b).
>
> I know that the x values are supposed to be x(1)=1000 and x(2)=0.5, as I previously solved this system by hand. However, fsolve gives me:
> x =
>
>  -73.3172
>    9.6259
>
> Is there anything I can do to improve he code so that it gives the correct answers, or something close?

When I pass the x values returned from fsolve back into f, I get some pretty small values:
>> x
x =
  -73.3172
    9.6259

>> f(x)
ans =
  -0.0000106970
   0.0000053239

When I pass your hand solution in, I get some much larger values:

>> f([1000; 0.5])
ans =
  -10.150
   39.951

So it looks like either the equations coded into your function don't match your theoretical equations or your hand solution is incorrect.

Tony Richardson



Reply | Threaded
Open this post in threaded view
|

RE: help fsolve not getting correct answers

Windhorn, Allen E [ACIM/LSA/MKT]
In reply to this post by Imane E
From: Help-octave [mailto:help-octave-bounces+allen.windhorn=[hidden email]] On Behalf Of Imane E

> I wrote this code:
> function y = f(x)
> ...
> And then used:
> [x, fval, info] = fsolve (@f, [1;2])
> to solve for x(1) and x(2) when a = b and c = d (where c and d are respective
> derivatives of a and b).

> I know that the x values are supposed to be x(1)=1000 and x(2)=0.5, as I
> previously solved this system by hand. However, fsolve gives me:
> x =
>  -73.3172
>   9.6259

If you make x = [1000;0.5] and put it into your function, the result is not zero:
>> x = [1000; 0.5];
>> f(x)
ans =
  -10.150
   39.951

So the function you have defined is not the same function you solved by
hand.  A solution to the Octave function is x = [-73.3172; 9.6259].
Probably the function you solved by hand is the correct one, since it comes
out to be round numbers, so you have to find out what is the difference
between that one and the one you wrote in Octave.  There are an infinite
number of functions that have the solution you are looking for, and we don't
have any way to guess which one you need, so you will have to troubleshoot
your code yourself.  This will be good practice.  Good luck to you.

P.S. The function may have more than one solution.

Regards,
Allen