# solving a very complicated equation implicitly

4 views (last 30 days)
Benjamin on 10 Apr 2019
Commented: Benjamin 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 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 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: Benjamin 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 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 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.