## solving a very complicated equation implicitly

### Benjamin (view profile)

on 10 Apr 2019
Latest activity Commented on by Benjamin

### Benjamin (view profile)

on 11 Apr 2019
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??

Benjamin

### Benjamin (view profile)

on 10 Apr 2019
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

### David Goodmanson (view profile)

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

### Benjamin (view profile)

on 11 Apr 2019
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?

### David Wilson (view profile)

Answer by David Wilson

### David Wilson (view profile)

on 11 Apr 2019
Edited by David Wilson

### David Wilson (view profile)

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:

Benjamin

### Benjamin (view profile)

on 11 Apr 2019
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)
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)

### David Wilson (view profile)

Answer by David Wilson

### David Wilson (view profile)

on 11 Apr 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.

Benjamin

### Benjamin (view profile)

on 11 Apr 2019
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.