Simultaneous Nonlinear Equations - Using fsolve - Specific Example (Troubleshoot)

Hi All,

I'm working on a simple function to solve a pair of nonlinear equations simultaneously using fsolve.

At the moment, I'm working on a specific example as a precursor to a more complex approach. Unfortunately, I seem to be having trouble getting the correct answer from my code. To save an aneurysm, I thought I'd ask for a troubleshoot.

The equations come from the Streeter-Phelps Model and are for the Critical Time & Critical Deficit. The equations are in the attached file in their original format.

My current Function looks like this:

function F=StreeterPhelps(k)

Kr=k(1);
Kd=k(2);

F(1)=((1/0.2775-Kr)*log(0.2775/Kr*(1-(1.173*(0.2775-Kr)/Kd*11))))-1.517;
F(2)=(Kd*11/0.2775*(0.2775/Kr*(1-(1.173*(0.2775-Kr)/Kd*11)))^-(Kr/0.2775-Kr))-6.643;

end

I'm trying to solve for Kr & Kd obviously but, at the moment, I'm getting Kr = 0.2894 and Kd = 0.2366 (Initial Guess [0.1;0.1]). The text I'm working from has the solution as 1.159 and 0.97 respectively. I'm wracking my brains but figure I have translated the equation badly?

Any help appreciated.

Answer 1

[John D'Errico]

You need to learn how to write an expression in MATLAB, in terms of the order of operations. If you don't ABSOLUTELY know how something will be parsed, then test it out! If you are not sure, then an extra pair of parens will not hurt, just make your code less readable.

For example, there is a difference between these two expressions:

1/0.2775-Kr

1/(0.2775-Kr)

I think you would agree with that claim? In the first case, MATLAB divides 1 by 0.2775, then it subtracts Kr. In the latter case, it subtracts Kr from 0.2775, then it divides that result into 1. Look at your code, then look the equations in that ODF. Which form do you really want to use?

You do similar things in several places in those equations.

[Samuel Dickinson]

Hi John,

Thanks for taking the time to reply. I certainly agree with everything you've said - writing the expressions into MatLab is not something I've had that much practice with.

I went through the function and rewrote it to:

F(1)=1/(0.2775-Kr)*log((0.2775/Kr)*(1-((1.173*(0.2775-Kr))/(Kd*11))))-1.517
F(2)=((Kd*11)/0.2775)*((0.2775/Kr)*(1-((1.173*(0.2775-Kr))/(Kd*11)))^-(Kr/0.2775-Kr))-6.643;

I'm getting Kr=1.1535 & Kd=0.9289. I think I'm getting close but perhaps still an issue or two somewhere?

Edit - A little more checking and I have it sorted now. Thanks again.