solving a very complicated equation implicitly
조회 수: 5 (최근 30일)
이전 댓글 표시
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213040/image.png)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213010/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213010/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213012/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213013/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213013/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213041/image.png)
Now once I show the function of g, it is obvious that
will need to be solved implicity.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213016/image.png)
Here is ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213012/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213012/image.png)
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213016/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213016/image.png)
This is the equation for
:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213013/image.png)
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
David Goodmanson
2019년 4월 11일
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.
답변 (2개)
David Wilson
2019년 4월 11일
편집: David Wilson
2019년 4월 11일
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:
![test.png](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213107/test.png)
David Wilson
2019년 4월 11일
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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213426/image.png)
then I re-arranged to
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/213427/image.png)
which I now solve for Xm.
참고 항목
카테고리
Help Center 및 File Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!