solving a very complicated equation implicitly
Show older comments
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:
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.
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.
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 Comments
Walter Roberson
on 10 Apr 2019
You are not going to be able to obtain a closed form formula. Although g( r) is a polynomial, it is being divided by r inside the integral. The integral itself is closed form, but involves log(x_m) . When you go to solve the integral == 1 portion for x_m then you need to take the root of an expression that involves rpf to various powers (polynomial-like) but also involves several exp() terms for each power of rpf .
Even if you substitute in a particular numeric rpf value, you end up with a sum of exponentials to try to find the root of. In the example numeric rpf value I tried, there turned out to be 5 real roots and 2 imaginary roots. With the exp() of the root structure, two of the roots turned out to be > 1 (as implied by the structure of the integral.)
Walter Roberson
on 10 Apr 2019
What is the potential range for rpf ?
An example rpf would help.
I hopped over to Maple to find solutions. With rpf = 2, there were two real solutions in the range you were looking for, namely 1.003781592, and 1,126335754. With rpf = 1/2, the real solution in the range you were looking for was 1.401500353
If rpf is constrained to positive, then it looks like there are solutions in the range of rpf slightly less than 1/4, to rpf = 1.. though I am still cross-checking that.
Walter Roberson
on 10 Apr 2019
Looks like the 1.126 etc. for rpf=2 was a false root due to insufficient precision in the calculations.
David Goodmanson
on 11 Apr 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.
Answers (2)
David Wilson
on 11 Apr 2019
Edited: David Wilson
on 11 Apr 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 Comment
David Wilson
on 11 Apr 2019
0 votes
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.
Categories
Find more on Numerical Integration and Differentiation in Help Center and File Exchange
on 11 Apr 2019
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
