Iteration of multiple nonlinear functions

4 次查看（过去 30 天）
Erik Eriksson
Erik Eriksson 2023-1-20
评论： Sam Chak 2024-5-30
I have an assignment where I need to clalculate the temperature distribution in a window with elctric heaters installed at nine different points.
I have made the correct temperature profile equation at each point. Now I need using iteration to find the temperature at a given time, and that is where my problem starts.
My code looks like this so far:
%Given values
k = 0.84;
Qh = 25;
Ti = -3;
To = -3;
hi = 6;
ho = 20;
dx = 0.002;
dy = 0.01;
n_steps = 1000;
tau = 0.078;
T1 = zeros(n_steps,1);
T1(1,1) = Ti;
T2 = zeros(n_steps,1);
T2(1,1) = Ti;
T3 = zeros(n_steps,1);
T3(1,1) = Ti;
T4 = zeros(n_steps,1);
T4(1,1) = Ti;
T5= zeros(n_steps,1);
T5(1,1) = Ti;
T6 = zeros(n_steps,1);
T6(1,1) = Ti;
T7 = zeros(n_steps,1);
T7(1,1) = Ti;
T8 = zeros(n_steps,1);
T8(1,1) = Ti;
T9 = zeros(n_steps,1);
T9(1,1) = Ti;
for i = 1:(n_steps-1)
syms T1 T2 T_ T4 T5 T6 T7 T8 T9
eqn1 = T1 == (1-2*hi*dy*tau/k-10.4*tau)*T1(i) + 2*hi*dy*tau*Ti/k + 10*tau*T2(i) + 0.4*tau*T4(i);
eqn2 = T2 == (1-10.4*tau)*T2(i) + 5*(T1(i)+T3(i))*tau + 0.4*tau*T5(i);
eqn3 = T3 == (1-2*ho*dy*tau/k-10.4*tau)*T3(i) + 2*ho*dy*tau*To/k + 10*tau*T2(i)+0.4*tau*T6(i);
eqn4 = T4 == (1-2*hi*dy*tau/k-10.4*tau)*T4(i) + 2*hi*dy*tau*Ti/k + 0.2*tau*(T1(i)+T7(i)) + 10*tau*T5(i);
eqn5 = T5 == (1-10.4*tau)*T5(i) + 0.2*(T2(i)+T8(i))*tau + 5*(T4(i)+T6(i))*tau;
eqn6 = T6 == (1-2*ho*dy*tau/k-10.4*tau)*T6(i) + 2*ho*dy*tau*To/k+0.2*(T3(i)+T9(i))*tau + 10*tau*T5(i);
eqn7 = T7 == (1-2*hi*dy*tau/k-10.4*tau)*T7(i) + 2*Qh*tau/k+2*hi*dy*tau*Ti/k+0.4*tau*T4(i) + 10*tau*T8(i);
eqn8 = T8 == (1-10.4*tau)*T8(i) + 0.4*tau*T5(i) + 5*(T7(i)+T9(i))*tau;
eqn9 = T9 == (1-2*ho*dy*tau/k-10.4*tau)*T9(i) + 2*ho*dy*tau*To/k + 0.4*tau*T6(i) + 10*T8(i)*tau;
s = solve([eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7 eqn8 eqn9],[T1 T2 T3 T4 T5 T6 T7 T8 T9])
s = struct2array(s);
s = double(s);
T1(i+1,1) = s(1);
T2(i+1,1) = s(2);
T3(i+1,1) = s(3);
T4(i+1,1) = s(4);
T5(i+1,1) = s(5);
T6(i+1,1) = s(6);
T7(i+1,1) = s(7);
T8(i+1,1) = s(8);
T9(i+1,1) = s(9);
end
%Temperature after 15 min
res_15min = [T1(225,1);T2(225,1);T3(225,1);T4(225,1);T5(225,1);T6(225,1);T7(225,1);T8(225,1);T9(225,1)]
%Temperature at steady-state
res_ss = [T1(525,1);T2(525,1);T3(525,1);T4(525,1);T5(525,1);T6(525,1);T7(525,1);T8(525,1);T9(525,1)]
I am unsure about the command "solve".
Does anyone know a better way of iterating these equation and then later plot them?

请先登录，再回答此问题。

回答（2 个）

Torsten
Torsten 2023-1-20
Ran in:
%Given values
k = 0.84;
Qh = 25;
Ti = -3;
To = -3;
hi = 6;
ho = 20;
dx = 0.002;
dy = 0.01;
n_steps = 1000;
tau = 0.078;
T1 = zeros(n_steps,1);
T1(1,1) = Ti;
T2 = zeros(n_steps,1);
T2(1,1) = Ti;
T3 = zeros(n_steps,1);
T3(1,1) = Ti;
T4 = zeros(n_steps,1);
T4(1,1) = Ti;
T5 = zeros(n_steps,1);
T5(1,1) = Ti;
T6 = zeros(n_steps,1);
T6(1,1) = Ti;
T7 = zeros(n_steps,1);
T7(1,1) = Ti;
T8 = zeros(n_steps,1);
T8(1,1) = Ti;
T9 = zeros(n_steps,1);
T9(1,1) = Ti;
for i = 1:(n_steps-1)
T1(i+1,1) = (1-2*hi*dy*tau/k-10.4*tau)*T1(i) + 2*hi*dy*tau*Ti/k + 10*tau*T2(i) + 0.4*tau*T4(i);
T2(i+1,1) = (1-10.4*tau)*T2(i) + 5*(T1(i)+T3(i))*tau + 0.4*tau*T5(i);
T3(i+1,1)= (1-2*ho*dy*tau/k-10.4*tau)*T3(i) + 2*ho*dy*tau*To/k + 10*tau*T2(i)+0.4*tau*T6(i);
T4(i+1,1)= (1-2*hi*dy*tau/k-10.4*tau)*T4(i) + 2*hi*dy*tau*Ti/k + 0.2*tau*(T1(i)+T7(i)) + 10*tau*T5(i);
T5(i+1,1)= (1-10.4*tau)*T5(i) + 0.2*(T2(i)+T8(i))*tau + 5*(T4(i)+T6(i))*tau;
T6(i+1,1)= (1-2*ho*dy*tau/k-10.4*tau)*T6(i) + 2*ho*dy*tau*To/k+0.2*(T3(i)+T9(i))*tau + 10*tau*T5(i);
T7(i+1,1)= (1-2*hi*dy*tau/k-10.4*tau)*T7(i) + 2*Qh*tau/k+2*hi*dy*tau*Ti/k+0.4*tau*T4(i) + 10*tau*T8(i);
T8(i+1,1)= (1-10.4*tau)*T8(i) + 0.4*tau*T5(i) + 5*(T7(i)+T9(i))*tau;
T9(i+1,1)= (1-2*ho*dy*tau/k-10.4*tau)*T9(i) + 2*ho*dy*tau*To/k + 0.4*tau*T6(i) + 10*T8(i)*tau;
end
plot((1:n_steps).',[T1,T2,T3,T4,T5,T6,T7,T8,T9])
  2 个评论
Erik Eriksson
Erik Eriksson 2023-1-24
Very simple and beatuful solution, exactly what i needed!
Thank you Sir!
Sam Chak
Sam Chak 2024-5-30
Hi @Erik Eriksson
If the beautiful solution you needed has been provided, please consider clicking 'Accept' ✔ on the answer to close the issue. Additionally, voting for @Alan Stevens' helpful answer serves as tokens of appreciation.

请先登录，再进行评论。

Alan Stevens
Alan Stevens 2023-1-20
Ran in:
Here's a quick and dirty way. I'll leave you to modify the following to record all the temperatures at every step.
%Given values
k = 0.84;
Qh = 25;
Ti = -3;
To = -3;
hi = 6;
ho = 20;
dx = 0.002;
dy = 0.01;
n_steps = 1000;
tau = 0.078;
V = [2*hi*dy*tau*Ti/k; 0; 2*ho*dy*tau*To/k; 2*hi*dy*tau*Ti/k; 0; 2*ho*dy*tau*To/k;
2*Qh*tau/k+2*hi*dy*tau*Ti/k; 0; 2*ho*dy*tau*To/k];
M = [(1-2*hi*dy*tau/k-10.4*tau), 10*tau, 0, 0.4*tau, 0, 0, 0, 0, 0;
5*tau, (1-10.4*tau), 5*tau, 0, 0.4*tau, 0,0,0,0;
0, 10*tau, (1-2*ho*dy*tau/k-10.4*tau), 0, 0, 0.4*tau, 0,0,0;
0.2*tau, 0, 0, (1-2*hi*dy*tau/k-10.4*tau), 10*tau, 0,0.2*tau,0,0;
0, 0.2*tau, 0, 5*tau, (1-10.4*tau), 5*tau, 0,0,0;
0, 0, 0.2*tau, 0, 10*tau, (1-2*ho*dy*tau/k-10.4*tau), 0, 0, 0.2*tau;
0, 0, 0, 0.4*tau, 0, 0, (1-2*hi*dy*tau/k-10.4*tau), 10*tau, 0;
0, 0, 0, 0, 0.4*tau, 0, 5*tau, (1-10.4*tau), 5*tau;
0, 0, 0, 0, 0, 0.4*tau, 0, 10*tau, (1-2*ho*dy*tau/k-10.4*tau)];
T = Ti*ones(9,1);
for i = 1:(n_steps-1)
T = M*T + V;
if i == 225
disp('Temps after 15mins')
disp(T')
figure
plot(1:9,T,'o-')
end
if i == 525
disp('Steady-state temps')
disp(T')
figure
plot(1:9,T,'o-')
end
end
Temps after 15mins
3.3086 3.2945 3.1006 6.0883 5.7174 5.6842 34.2074 29.9183 27.5795
Steady-state temps
3.6452 3.6307 3.4263 6.3811 6.0066 5.9675 34.5453 30.2559 27.9065

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by