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)
采纳的回答
Andrew Newell
2011-2-5
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?
5 个评论
hashkul
2011-2-5
Andrew Newell
2011-2-5
Was your comment before or after my edit?
hashkul
2011-2-5
hashkul
2011-2-5
Andrew Newell
2011-2-5
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?
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Calculus 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
