How to solve non-linear equations with integrations?

I want to solve a group of equations like this:

f1 = int(pl*x,x,-L,0)+int(pr*x,x,0,L)+M/d;

f2 = int(pl,x,-L,0)*sin(alph+tilt)+int(pr,x,0,L)*sin(alph-tilt)-G/d-2*L*P_e*sin(alph)*cos(tilt);

f3 = int(pl,x,-L,0)*cos(alph+tilt)-int(pr,x,0,L)*cos(alph-tilt)+2*L*P_e*sin(alph)*sin(tilt);

in which the pl and pr are:

pl = -12*mu_L*rho_S*(int((Uxl*x/(deltal^3)),x,0,x)+Cl*int((1/(deltal^3)),x,0,x))/rho_L+P0;

pr = -12*mu_L*rho_S*(int((Uxr*x/(deltar^3)),x,0,x)+Cr*int((1/(deltar^3)),x,0,x))/rho_L+P0;

The solving prt of code like:

initial_guess = [0.5, 0.1, 2];

solve_equations = @(vars) double(subs([f1, f2, f3], {x1,y1,Uratio}, vars));

result = fsolve(solve_equations, initial_guess);

x1 = result(1);

y1 = result(2);

Uratio = result(3);

But the calculating process time was extremly long and couldn't get the solution, I thought the reason is: there are lots of integrations in the equations and the fsolve function cannot solve the int() part. Is this correct? How to solve it?

[Torsten]

Where are x1,y1 and Uratio in your three equations ? And the equations to solve should be f1 = 0, f2 = 0 and f3 = 0 ?

[Sam Chak]

Hi @Yuting

Are the equations related to Mass Flow Rate and Heat Transfer?

Do the variables alph and tilt denote α and θ, respectively?

[Yuting]

Hi @Sam Chak,

yes the euqations represent to a close-contact melting problem. In which the alpha donates α as the a predefined Structural parameter, and tilt donates the inclination angle θ of the heat source

[Yuting]

Hi @Torsten,

Where are x1,y1 and Uratio in your three equations ?

The x1,y1 and Uratio are included in the expression of Uxl, Uxr, Cl, Cr, deltal, deltar, and P0.

And the equations to solve should be f1 = 0, f2 = 0 and f3 = 0 ?

Yes~ Excatly~

[Torsten]

Do some of your "parameters" Uxl, Uxr, Cl, Cr, deltal, deltar, P0, mu_L, rho_S, alph, tilt, P_e ... also depend on x ?

[Yuting]

@Torsten No, the other parameters like mu_L, rho_S, alph, tilt, P_e are not related to x~

[Torsten]

Then follow @Star Strider to work with "matlabFunction", but compute pl and pr first since the result will also depend on x. This will influence int(pl*x,x,-L,0) etc. in the computation of f1, f2 and f3.

syms x L M d alph tilt G P_e delta rho_L Uxl Uxr rho_S mu_L P0 Cr Cl deltal deltar 
pl = -12*mu_L*rho_S*(int((Uxl*x/(deltal^3)),x,0,x)+Cl*int((1/(deltal^3)),x,0,x))/rho_L+P0;
pr = -12*mu_L*rho_S*(int((Uxr*x/(deltar^3)),x,0,x)+Cr*int((1/(deltar^3)),x,0,x))/rho_L+P0;
f1 = int(pl*x,x,-L,0)+int(pr*x,x,0,L)+M/d

f1 = 

f2 = int(pl,x,-L,0)*sin(alph+tilt)+int(pr,x,0,L)*sin(alph-tilt)-G/d-2*L*P_e*sin(alph)*cos(tilt)

f2 = 

f3 = int(pl,x,-L,0)*cos(alph+tilt)-int(pr,x,0,L)*cos(alph-tilt)+2*L*P_e*sin(alph)*sin(tilt)

f3 = 

equations = matlabFunction(f1,f2,f3)

equations = function_handle with value:
    @(Cl,Cr,G,L,M,P0,P_e,Uxl,Uxr,alph,d,deltal,deltar,mu_L,rho_L,rho_S,tilt)deal(M./d+(L.^2.*1.0./deltal.^3.*(L.^2.*Uxl.*mu_L.*rho_S.*3.0-Cl.*L.*mu_L.*rho_S.*8.0))./(rho_L.*2.0)-(L.^2.*1.0./deltar.^3.*(L.^2.*Uxr.*mu_L.*rho_S.*3.0+Cr.*L.*mu_L.*rho_S.*8.0))./(rho_L.*2.0),sin(alph+tilt).*(L.*P0-(L.*1.0./deltal.^3.*(L.^2.*Uxl.*mu_L.*rho_S.*2.0-Cl.*L.*mu_L.*rho_S.*6.0))./rho_L)-G./d+sin(alph-tilt).*(L.*P0-(L.*1.0./deltar.^3.*(L.^2.*Uxr.*mu_L.*rho_S.*2.0+Cr.*L.*mu_L.*rho_S.*6.0))./rho_L)-L.*P_e.*sin(alph).*cos(tilt).*2.0,cos(alph+tilt).*(L.*P0-(L.*1.0./deltal.^3.*(L.^2.*Uxl.*mu_L.*rho_S.*2.0-Cl.*L.*mu_L.*rho_S.*6.0))./rho_L)-cos(alph-tilt).*(L.*P0-(L.*1.0./deltar.^3.*(L.^2.*Uxr.*mu_L.*rho_S.*2.0+Cr.*L.*mu_L.*rho_S.*6.0))./rho_L)+L.*P_e.*sin(alph).*sin(tilt).*2.0)

[Yuting]

@Torsten Thank you so much!!! Very important!!

Answer 1

[Star Strider]

It is not possible to solve this numerically without substituting numerical values for the known variables.

That aside, use the matlabFunction function to create an anonymous function that fsolve can work with —

syms pl x pr L M d alph tilt G P_e delta rho_L Uxl Uxr rho_S mu_L P0 Cr Cl deltal deltar 

f1 = int(pl*x,x,-L,0)+int(pr*x,x,0,L)+M/d

f1 = 

f2 = int(pl,x,-L,0)*sin(alph+tilt)+int(pr,x,0,L)*sin(alph-tilt)-G/d-2*L*P_e*sin(alph)*cos(tilt)

f2 = 

f3 = int(pl,x,-L,0)*cos(alph+tilt)-int(pr,x,0,L)*cos(alph-tilt)+2*L*P_e*sin(alph)*sin(tilt)

f3 = 

pls = -12*mu_L*rho_S*(int((Uxl*x/(deltal^3)),x,0,x)+Cl*int((1/(deltal^3)),x,0,x))/rho_L+P0;
prs = -12*mu_L*rho_S*(int((Uxr*x/(deltar^3)),x,0,x)+Cr*int((1/(deltar^3)),x,0,x))/rho_L+P0;

f1 = simplify(subs(f1, {pl,pr},{pls,prs}), 500)

f1 = 

f2 = simplify(subs(f1, {pl,pr},{pls,prs}), 500)

f2 = 

f3 = simplify(subs(f1, {pl,pr},{pls,prs}), 500)

f3 = 

equations = matlabFunction(f1,f2,f3)

equations = function_handle with value:
    @(Cl,Cr,L,M,Uxl,Uxr,d,deltal,deltar,mu_L,rho_L,rho_S,x)deal(M./d+(L.^2.*1.0./deltal.^3.*1.0./deltar.^3.*mu_L.*rho_S.*x.*(Cl.*deltar.^3.*2.0-Cr.*deltal.^3.*2.0+Uxl.*deltar.^3.*x-Uxr.*deltal.^3.*x).*3.0)./rho_L,M./d+(L.^2.*1.0./deltal.^3.*1.0./deltar.^3.*mu_L.*rho_S.*x.*(Cl.*deltar.^3.*2.0-Cr.*deltal.^3.*2.0+Uxl.*deltar.^3.*x-Uxr.*deltal.^3.*x).*3.0)./rho_L,M./d+(L.^2.*1.0./deltal.^3.*1.0./deltar.^3.*mu_L.*rho_S.*x.*(Cl.*deltar.^3.*2.0-Cr.*deltal.^3.*2.0+Uxl.*deltar.^3.*x-Uxr.*deltal.^3.*x).*3.0)./rho_L)

I leave the rest to you.

.

[Yuting]

Thanks a lot !!!! I will try this way!

[Yuting]

Hello~ @Star Strider,

I have a question about the code f1 = simplify(subs(f1, {pl,pr},{pls,prs}), 500),

what does the "500" mean?

[Star Strider]

My pleasure!

The ‘500’ tells simplify too keep simplifying until it cannot simplify the expressins further, or reaches the limit of 500 iterations. Usually 500 is enough, however if not specifically stated, simplify stops after one iteration. For complicated expressions, that is rarely enough.