Iteration of multiple nonlinear functions

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?

Answer 1

[Torsten]

%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])

[Erik Eriksson]

Very simple and beatuful solution, exactly what i needed!

Thank you Sir!

[Sam Chak]

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.

Answer 2

[Alan Stevens]

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