Using QUAD function in my code

Dear all,

I wanted to ask if my implementation of the quad function is right in this code I have laid out? I am getting negative answers (negative volumes) and they most definitely need to be positive.

However I have no idea where I am going wrong?

The idea is to find V from the following differential equation:

dF3MEK/dV = r1-2*r2

So when we re-arrange we get:

dF3MEK/(r1-2*r2) = V

Therefore in order to integrate with respect to x on the L.H.S I need to determine the following integral (as also seen in the code at line 20).

dF3MEK/dx = differential

Everything is in symbolic x form till now and then I convert the resulting equation using matlabFunction into the correct array formatting required for quad.

And then I run quad on the equation to determine the integral for which I should get a positive volume. I however do not know what the problem is at the moment and why I am getting a negative volume.

Thank you very much in advance. The code is attached below:

          x = sym('x', 'real');
          Ptot = 200; %kPa
          T = 450; %K
          x1H2O = 0.102398524;
          x1 = 0.22;  
          F2SBA = 107.5882086; %mol/s
          S1 = -3927667667.6679700000*x^6 + 7615096701.2234400000*x^5 - 5984633718.0112800000*x^4 + 2425863897.0966400000*x^3 - 529294313.9777970000*x^2 + 57463477.2858973000*x - 2227467.2236688300;
          F3SBA = F2SBA*(1-x);
          F3H2 = F2SBA*(x);
          F3MEK = (F2SBA*x*S1)/(S1+2)
          differential = diff(F3MEK);
          F3BP = (F2SBA*x)/(S1+2)
          F3H2O = (F2SBA*x)/(S1+2)+(F2SBA/(1-x1H2O))*x1H2O
          F3 = F3SBA + F3H2 + F3MEK + F3H2O + F3BP
          PSBA = ((F3SBA)/(F3))*Ptot;
          PH2 = ((F3H2)/(F3))*Ptot;
          PMEK = ((F3MEK)/(F3))*Ptot;
          rxn1 = ((8.290*10^5)*exp(-6403/T)*(PSBA-((PMEK*PH2)/((3.538*10^8)*exp(-7100/T)))))/((1+(8.804*10^-5)*exp(3298/T)*PMEK)^2)
          rxn2 = (1.22*10^2)*exp(-9860/T)*(PMEK^2)
          integrand = (differential)/(rxn1-2*rxn2);
          quadfunction = matlabFunction(integrand)
          f = @(x)quadfunction;
          Volume = quad(f, 0, 0.22)

Answer 1

[Andrew Newell]

It looks like you have a singularity in the integrand at x=0.08. Try this command to look at it:

ezplot(integrand,[0 0.22])

If you can determine the exact location of the singularity, you could try integrating the intervals below and above using quadgk.

EDIT: I looked at this expression in Maple, and there is indeed a singularity at 0.07965636277, as well as about 30 other complex values. My suggestion above didn't work. I think your real problem is that you have a pretty nasty rational function with powers of x up to 28 in the denominator and coefficients spanning about 40 orders of magnitude. What you really need to do is think carefully about this problem - can it be made simpler?

[hashkul]

That is true however it doesn't seem to be making a difference even if I integrate till 0.07. It always brings a negative output, and from the plot it is clear it is not a negative area under the curve. So I actually still do not understand what is happening?

[Andrew Newell]

Was your comment before or after my edit?

[hashkul]

Thanks for your prompt reply! My comment was before your edit, and it does make sense that the might have been made a bit too complex. I will have to look into this and try to cut it down I guess.

I believe the equations are as they are, and they can't really be cut down. I am failing to see how this could be made simple really.

Is there any way to use trapz method with this?

[hashkul]

Second sentence is meant to read: I believe the equations are going to remain as they are...

[Andrew Newell]

Here's one idea. Aside from the singularity, the curve looks quite smooth. You could sample some points, avoiding the area near 0.08, fit them with a polynomial and then integrate that. It depends how accurate your answer needs to be. But also ask yourself - does the singularity make sense in terms of your problem?