solving a very complicated equation implicitly

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

[Walter Roberson]

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.)

[Benjamin]

@WalterRoberson My goal is just to find the root that is between the values of 1 and 2, so it could be bounded like that. Were any of the roots you found between those values? If so, then I believe you have exactly what I need.

Is there any chance you could post the code you tried so that I could play around with it? Also, I took the integral of the function times 1/r in mathematica and simplfied it and got this below. I know that it is hard to read but basically everything is separated out into powers of x, which can be seen at the end of each block.

Also, could you explain where the exp's are coming from exactly? I think I am confused by what you mean. I just want to use MATLAB to find the roots of the equation, which I think you did. Could you please post your code so I can see what you are doing?

1/60 (-60 C10 + 90 C20 - 110 C30 + 125 C40 - 137 C50 +

rpf (-60 C11 + 90 C21 - 110 C31 + 125 C41 - 137 C51 +

rpf (-60 C12 + 90 C22 - 110 C32 + 125 C42 - 137 C52 +

rpf (-60 C13 + 90 C23 - 110 C33 + 125 C43 - 137 C53 -

rpf (60 C14 - 90 C24 + 110 C34 - 125 C44 +

137 C54 + (60 C15 - 90 C25 + 110 C35 - 125 C45 +

137 C55) rpf)))) +

60 (C10 - 2 C20 + 3 C30 - 4 C40 + 5 C50 +

rpf (C11 - 2 C21 + 3 C31 - 4 C41 + 5 C51 +

rpf (C12 - 2 C22 + 3 C32 - 4 C42 + 5 C52 +

rpf (C13 - 2 C23 + 3 C33 - 4 C43 + 5 C53 +

rpf (C14 - 2 C24 + 3 C34 - 4 C44 +

5 C54 + (C15 - 2 C25 + 3 C35 - 4 C45 +

5 C55) rpf))))) x +

30 (C20 - 3 C30 + 6 C40 - 10 C50 +

rpf (C21 - 3 C31 + 6 C41 - 10 C51 +

rpf (C22 - 3 C32 + 6 C42 - 10 C52 +

rpf (C23 - 3 C33 + 6 C43 - 10 C53 +

rpf (C24 - 3 C34 + 6 C44 -

10 C54 + (C25 - 3 C35 + 6 C45 -

10 C55) rpf))))) x^2 +

20 (C30 - 4 C40 + 10 C50 +

rpf (C31 - 4 C41 + 10 C51 +

rpf (C32 - 4 C42 + 10 C52 +

rpf (C33 - 4 C43 + 10 C53 +

rpf (C34 - 4 C44 +

10 C54 + (C35 - 4 C45 + 10 C55) rpf))))) x^3 +

15 (C40 - 5 C50 +

rpf (C41 - 5 C51 +

rpf (C42 - 5 C52 +

rpf (C43 - 5 C53 +

rpf (C44 - 5 C54 + C45 rpf - 5 C55 rpf))))) x^4 +

12 (C50 +

rpf (C51 + rpf (C52 + rpf (C53 + rpf (C54 + C55 rpf))))) x^5 +

60 (C00 - C10 + C20 - C30 + C40 - C50 +

rpf (C01 - C11 + C21 - C31 + C41 - C51 +

rpf (C02 - C12 + C22 - C32 + C42 - C52 +

rpf (C03 - C13 + C23 - C33 + C43 - C53 +

rpf (C04 - C14 + C24 - C34 + C44 -

C54 + (C05 - C15 + C25 - C35 + C45 -

C55) rpf))))) Log[x])

[Walter Roberson]

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]

Looks like the 1.126 etc. for rpf=2 was a false root due to insufficient precision in the calculations.

[Benjamin]

@WalterRoberson The value of rpf has to be between 0 and 1, and Xm which is solved for needs to be between 1 and 2.

Ideally, I'd like to vectorize it and numerically solve for Xm for each value of rpf. So a good place to start would be to test rpf = 0, 0.1, 0.2 ... 1.

[Benjamin]

Just to follow up on the last comment as I see you posted at the same time, I only care about rpf's between 0 and 1 and Xm values between 1 and 2.

[Benjamin]

Sorry, I just realized a mistake. Outside of the integral there should be a sqrt(2)*pi/3. I have fixed the original question so you can see this mistake. I'm really sorry about that and appreciate your help.Basically, I want to get a lot of Xm values for different rpf values, and fit that curve.

Edit: @WalterRoberson Also, I think for an rpf = 0.5, then Xm should approximately be between 1.55 and 1.6. I could be wrong about that, but if that is what you get with accounting for the sqrt(2)*pi/3, then it must be correct. Do you think MATLAB can help me with this, or should I download Maple?

[David Goodmanson]

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.

[Benjamin]

I had to update the g_0 function (see above). Thanks for your help here. I can visually see where it crosses, can that be automated to return the actual number? I'm not sure how to vectorize this code for many values of rpf. Note that I had to edit the coefficients. The new coefficients have been updated in the original question.

Also, could you post your code with cubic splines as an answer? It would be great if it could be vectorized and just outputs where it crosses. I checked your code again with my new coeffs and everything, and I get the answer I expect! I would very much appreciate it if you could post it as an answer so I can accept it. Also, currenlty it will take a long time to check where each rpf crosses 1. How to automate?

Answer 1

[David Wilson]

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:

[Benjamin]

Wow thanks for this answer, i'm still working through what you did. Unfortunately, I need to change my g_0 function from:

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; 

to substract out the 1, divide by (4*rpf) and divide by (1/1.55)

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); 

However, when I do this, it just keeps throwing an error. If you change it in your code, does it run for you? I cannot get it to run. Your code works great, I guess just not for my corrected g_0 function, I can't figure out what's going wrong. Note that I had to update the coefficients since the function changed. They are now corrected in the original question.

Also why do you subtract out 3/sqrt(2)/pi? I can't seem to follow the algebra that would cause me to substract it. Looking at the integral, shouldn't it be divided out?

It currently gives me this:

Error in fsolve (line 408)

        trustnleqn(funfcn,x,verbosity,gradflag,options,defaultopt,f,JAC,...

Error in xm_solver_2>@(rpf)fsolve(@(Xm)Q(Xm,rpf)-3/(sqrt(2)*pi*g_0(rpf)),Xm0)

Error in xm_solver_2 (line 77)
Xm = arrayfun(Xm,rpf)

Answer 2

[David Wilson]

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.

[Benjamin]

I tried dot dividing first and that did not work. The error gets thrown when I remove the 1 from the equation. If you put it in your code, you can see that the error is thrown with or without dot dividing.