22 views (last 30 days)

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:

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

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:

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.

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!