Solve taking too long, how to optimize code?

6 visualizaciones (últimos 30 días)
Brock Ferrari
Brock Ferrari el 26 de Sept. de 2021
Comentada: Walter Roberson el 27 de Sept. de 2021
Basically, what the title asks. Code is here:
k = 180;
length = .05;
b = .01;
T0 = 180;
Tinf = 25;
h = 25;
Tsurr = 290;
epsilon = 0.9;
theta = atan((b/2)/length);
stef = 5.67*10^(-8);
xchange = b;
syms T1 T2 T3 T4 T5
eqn1 = (1-0.5*(xchange/length))*(T0 - T1) + (1-1.5*(xchange/length))*(T2-T1) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T1) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T1+273)^4) == 0;
eqn2 = (1-1.5*(xchange/length))*(T1 - T2) + (1-2.5*(xchange/length))*(T3-T2) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T2) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T2+273)^4) == 0;
eqn3 = (1-2.5*(xchange/length))*(T2 - T3) + (1-3.5*(xchange/length))*(T4-T3) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T3) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T3+273)^4) == 0;
eqn4 = (1-3.5*(xchange/length))*(T3 - T4) + (1-4.5*(xchange/length))*(T5-T4) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T4) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T4+273)^4) == 0;
eqn5 = 2*k*(xchange/2)*tan(theta)*((T4-T5)/xchange) + 2*h*((xchange/2)/cos(theta))*(Tinf - T5) + 2*epsilon*stef*((xchange/2)/cos(theta))*(Tsurr^4 - (T5 + 273)^4 ) == 0;
sol = solve([eqn1, eqn2, eqn3, eqn4, eqn5], [T1, T2, T3, T4, T5]);
T1sol = double(sol.T1)
T2sol = double(sol.T2)
T3sol = double(sol.T3)
T4sol = double(sol.T4)
T5sol = double(sol.T5)
Also, how long would this theoretically take to solve?
Thanks.

Respuestas (3)

Matt J
Matt J el 26 de Sept. de 2021
Well, one solution is clearly,
T1=T2=T3=T4=T5=(Tsurr-273)
as you can verify by direct subsitution.
  1 comentario
Walter Roberson
Walter Roberson el 27 de Sept. de 2021
This is not one of the 1024 roots that my code comes up with.

Iniciar sesión para comentar.


Alan Weiss
Alan Weiss el 27 de Sept. de 2021
Editada: Alan Weiss el 27 de Sept. de 2021
Another point of view in case Matt's solution is not satisfactory: work completely in floating point computation rather than symbolic math computation. If you use the Problem-Based Optimization Workflow, you will basically have the same equations as before, simply using optimization variables rather than symbolic variables.
Good luck,
Alan Weiss
MATLAB mathematical toolbox documentation

Walter Roberson
Walter Roberson el 27 de Sept. de 2021
Editada: Walter Roberson el 27 de Sept. de 2021
k = 180;
length = .05;
b = .01;
T0 = 180;
Tinf = 25;
h = 25;
Tsurr = 290;
epsilon = 0.9;
theta = atan((b/2)/length);
stef = 5.67*10^(-8);
xchange = b;
syms T1 T2 T3 T4 T5
eqn1 = (1-0.5*(xchange/length))*(T0 - T1) + (1-1.5*(xchange/length))*(T2-T1) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T1) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T1+273)^4) == 0;
eqn2 = (1-1.5*(xchange/length))*(T1 - T2) + (1-2.5*(xchange/length))*(T3-T2) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T2) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T2+273)^4) == 0;
eqn3 = (1-2.5*(xchange/length))*(T2 - T3) + (1-3.5*(xchange/length))*(T4-T3) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T3) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T3+273)^4) == 0;
eqn4 = (1-3.5*(xchange/length))*(T3 - T4) + (1-4.5*(xchange/length))*(T5-T4) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T4) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T4+273)^4) == 0;
eqn5 = 2*k*(xchange/2)*tan(theta)*((T4-T5)/xchange) + 2*h*((xchange/2)/cos(theta))*(Tinf - T5) + 2*epsilon*stef*((xchange/2)/cos(theta))*(Tsurr^4 - (T5 + 273)^4 ) == 0;
E1 = [eqn1, eqn2, eqn3, eqn4, eqn5];
%{
partial2 = solve(E1(1), T2)
E2 = subs(E1(2:end), T2, partial2)
partial3 = solve(E2(1), T3)
E3 = subs(E2(2:end), T3, partial3)
partial4 = solve(E3(1), T4)
E4 = subs(E3(2:end), T4, partial4)
partial5 = solve(E4(1), T5)
E5 = subs(E4(2:end), T5, partial5)
%}
partials = solve(E1(1:4), [T2, T3, T4, T5]);
E5 = subs(lhs(E1(end))-rhs(E1(end)), partials);
At this point you can expand(E5) to get a plain polynomial in T1.
According to Maple, that plain polynomial is in degree 1024.
You should be able to take the following steps, but they take too much time for me to be able to run them and show the output here
tic
P = sym2poly(E5);
R = roots(P);
mask = imag(R) == 0;
T1_sol = R(mask);
T2_sol = double(subs(partials.T2, T1, T1_sol));
T3_sol = double(subs(partials.T3, T1, T1_sol));
T4_sol = double(subs(partials.T4, T1, T1_sol));
T5_sol = double(subs(partials.T5, T1, T1_sol));
toc
T1_sol
T2_sol
T3_sol
T4_sol
T5_sol
  1 comentario
Walter Roberson
Walter Roberson el 27 de Sept. de 2021
T1_sol =
-9069.32334661296
-919.882898288149
177.024444908286
-24.4166503435627
T2_sol =
27738.4267958556
-2336.41842690512
174.081315968202
-287.461625512644
T3_sol =
7095829.16820471
-4126.23669567827
171.16801338191
-657.54975176523
T4_sol =
4.81617135399466e+16
-2960.81576251875
168.277608495836
-1280.43369676003
T5_sol =
3.06585398161644e+56
3425.69627315686
165.36383335676
-3127.23531356248

Iniciar sesión para comentar.

Categorías

Más información sobre Programming en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by