MATLAB: Using QUAD function in the code

differential equationsdifferentiationintegrationMATLABmatlabfunctionquadsymsymbolicSymbolic Math Toolboxsyms

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)

Best Answer

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?
Related Question