solving a very complicated equation implicitly
7 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hello, I have a very complicated equation that needs to be integrated and solved implicitly. I am not sure if MATLAB can handle this, so I thought I would ask here for help.
I have this equation:
However, is somewhat complicated, but at the end of the day, it's just a polynomial. is actually comprised of two functions: and . However, has no dependence on r so it can be pulled out. Thus I have:
Now once I show the function of g, it is obvious that will need to be solved implicity.
Here is
g = C00*(r-1).^0*rpf^0 + C10*(r-1).^1*rpf^0 + C20*(r-1).^2*rpf^0 + C30*(r-1).^3*rpf^0 + ...
+ C40*(r-1).^4*rpf^0 + C50*(r-1).^5*rpf^0 + C01*(r-1).^0*rpf^1 + C11*(r-1).^1*rpf^1 + ...
+ C21*(r-1).^2*rpf^1 + C31*(r-1).^3*rpf^1 + C41*(r-1).^4*rpf^1 + C51*(r-1).^5*rpf^1 + ...
+ C02*(r-1).^0*rpf^2 + C12*(r-1).^1*rpf^2 + C22*(r-1).^2*rpf^2 + C32*(r-1).^3*rpf^2 + ...
+ C42*(r-1).^4*rpf^2 + C52*(r-1).^5*rpf^2 + C03*(r-1).^0*rpf^3 + C13*(r-1).^1*rpf^3 + ...
+ C23*(r-1).^2*rpf^3 + C33*(r-1).^3*rpf^3 + C43*(r-1).^4*rpf^3 + C53*(r-1).^5*rpf^3 + ...
+ C04*(r-1).^0*rpf^4 + C14*(r-1).^1*rpf^4 + C24*(r-1).^2*rpf^4 + C34*(r-1).^3*rpf^4 + ...
+ C44*(r-1).^4*rpf^4 + C54*(r-1).^5*rpf^4 + C05*(r-1).^0*rpf^5 + C15*(r-1).^1*rpf^5 + ...
+ C25*(r-1).^2*rpf^5 + C35*(r-1).^3*rpf^5 + C45*(r-1).^4*rpf^5 + C55*(r-1).^5*rpf^5;
Now this equation needs to be multiplied by 1/r and integrated with the corresponding limits. I would then like to solve for as a function of the rpf. Note that rpf goes from 0 to 0.7. So I want to solve implicitly for for different values of rpf; namely, rpf = 0:0.01:0.7.
This is the equation for :
g_0 = @(rpf) (3*rpf./(1-rpf)+(1)*A1*rpf.^1 + (2)*A2*rpf.^2 + (3)*A3*rpf.^3 + ...
(4)*A4*rpf.^4 + (10)*A10*rpf.^10 + (14)*A14*rpf.^14)/(4.*rpf)/(1/1.55);
and below are all the coefficients:
C12 = -7.057351859 ;
C13 = 11.25652658 ;
C14 = -27.94282684 ;
C15 = 7.663030197 ;
C22 = 30.19116817 ;
C23 = -110.0869103 ;
C24 = 195.580299 ;
C25 = 23.52636596 ;
C32 = -71.65721882 ;
C33 = 324.0458389 ;
C34 = -311.2587989 ;
C35 = -429.8278926 ;
C42 = 82.67213601 ;
C43 = -325.747121 ;
C44 = -127.5075666 ;
C45 = 1184.468297 ;
C52 = -30.4648185 ;
C53 = 52.49640407 ;
C54 = 467.0247693 ;
C55 = -1014.5236 ;
%Known Coefficients
C00 = 1 ;
C10 = 0 ;
C20 = 0 ;
C30 = 0 ;
C40 = 0 ;
C50 = 0 ;
C01 = 0 ;
C11 = -2.903225806 ;
C21 = 0.967741935 ;
C31 = 0.322580645 ;
C41 = 0 ;
C51 = 0 ;
C02 = 0 ;
C03 = 0 ;
C04 = 0 ;
C05 = 0 ;
A1 = -0.419354838709677;
A2 = 0.581165452653486;
A3 = 0.643876195271502;
A4 = 0.657898291526944;
A10 = 0.651980607965746;
A14 = -2.6497739490251;
Can anyone help me input this into MATLAB in a way that solves for X_M for the different values of rpf??
9 commentaires
David Goodmanson
le 11 Avr 2019
Hi Benjamin/Walter,
I decided to try my favorite method of plotting the integrand and integral so you can see what is going on.
constants = (what you have)
rpf = .5;
% upper limit on r below allows the indefinite integral
% to always be > 1 for 0 <= rpf <= .999
r = linspace(1,2-rpf,1000)
g_0 = (what you have)
g = (what you have)
f = (sqrt(2)*pi/3)*g_0*g./r;
figure(1)
plot(r,f)
grid on
Intf = cumtrapz(r,f)
figure(2)
plot(r,Intf)
grid on
For rpf = .5 the indefinite integral crosses 1 at approximately r = 1.226 = Xm. Despite the fact that cumtrapz is pretty crude, it's not bad here because the integrand is smooth. A better integration and interpolation method using cubic splines comes up with 1.2266.
With interpolation using Intf as the independent variable and r as the dependent variable, the cumtrapz approach could be automated to find Xm for many values of rpf.
Réponses (2)
David Wilson
le 11 Avr 2019
Modifié(e) : David Wilson
le 11 Avr 2019
OK, try this:
My approach was to embed an integral inside a root solver. I basically used your code and coefficients, but changed them to anonymous functions.
I needed a starting guess, and actually it worked for all your cases. I wrapped the entire thing up also to use arrayfun in an elegant way to do the whole rpf vector in a single line.
g_0 = @(rpf) 1 + 3*rpf./(1-rpf)+(1)*A1*rpf.^1 + (2)*A2*rpf.^2 + (3)*A3*rpf.^3 + ...
(4)*A4*rpf.^4 + (10)*A10*rpf.^10 + (14)*A14*rpf.^14;
g = @(r,rpf) C00*(r-1).^0*rpf^0 + C10*(r-1).^1*rpf^0 + C20*(r-1).^2*rpf^0 + C30*(r-1).^3*rpf^0 + ...
+ C40*(r-1).^4*rpf^0 + C50*(r-1).^5*rpf^0 + C01*(r-1).^0*rpf^1 + C11*(r-1).^1*rpf^1 + ...
+ C21*(r-1).^2*rpf^1 + C31*(r-1).^3*rpf^1 + C41*(r-1).^4*rpf^1 + C51*(r-1).^5*rpf^1 + ...
+ C02*(r-1).^0*rpf^2 + C12*(r-1).^1*rpf^2 + C22*(r-1).^2*rpf^2 + C32*(r-1).^3*rpf^2 + ...
+ C42*(r-1).^4*rpf^2 + C52*(r-1).^5*rpf^2 + C03*(r-1).^0*rpf^3 + C13*(r-1).^1*rpf^3 + ...
+ C23*(r-1).^2*rpf^3 + C33*(r-1).^3*rpf^3 + C43*(r-1).^4*rpf^3 + C53*(r-1).^5*rpf^3 + ...
+ C04*(r-1).^0*rpf^4 + C14*(r-1).^1*rpf^4 + C24*(r-1).^2*rpf^4 + C34*(r-1).^3*rpf^4 + ...
+ C44*(r-1).^4*rpf^4 + C54*(r-1).^5*rpf^4 + C05*(r-1).^0*rpf^5 + C15*(r-1).^1*rpf^5 + ...
+ C25*(r-1).^2*rpf^5 + C35*(r-1).^3*rpf^5 + C45*(r-1).^4*rpf^5 + C55*(r-1).^5*rpf^5;
Xm0 = 1;
Q = @(Xm,rpf) integral(@(r) g(r,rpf)./r, 1, Xm)
Xm = @(rpf) fsolve(@(Xm) Q(Xm, rpf)-3/(sqrt(2)*pi*g_0(rpf)), Xm0);
rpf = [0:0.01:0.7]';
Xm = arrayfun(Xm,rpf)
plot(rpf, Xm,'.-')
And the plot looks like:
1 commentaire
David Wilson
le 11 Avr 2019
A quick answer is you probably need to dot-divide (./) the term with the (vector) rpf.
2nd point. Because I used a root finder, which searches for f(x) = 0, then I had to convert your original task which was to find the value of Xm that made your integral = to 1, to an expression that equalled zero.
So I started with (in simplied nomenclature)
then I re-arranged to
which I now solve for Xm.
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!