Solve the differential equation using Crank-Nicolson method and Newon iterative method.

I'm tring to solve the eqation something like this to get the transient temperature of each node:

where, 𝑡 is time; 𝑃(𝑡) is the input power; 𝑁𝑖 is the note 𝑖 (𝑖 = 1,2,.. .,n+1) (𝑁1 is the heat source and 𝑁𝑛+1 is the constant temperature environment); 𝑇𝑖 = 𝑇𝑖(𝑡) is the temperature at 𝑁𝑖 while t; 𝑅𝑡ℎ𝑖 is the thermal resistance between 𝑁𝑖 and 𝑁𝑖+1; 𝐶𝑡ℎ 𝑖 is the thermal capacity of 𝑁𝑖.

With the initial condition:

I try to discretize in time with Crank Nicolson method and solve it using Newton iterative method:

But I don't know how to do that with f(t) equation contains Ti-1, Ti and Ti+1 together, could anyone know how to solve this? I would like to obtain Ti(t) , thanks a lot.

I have tried it like this, but it seems like a wrong code.

n = 4; 
dt = 0.01; 
timesteps = 1000; 
R = [6.9364, 8.3004, 2.8952, 12.0162]; 
C = [1, 2, 1, 2]; 
P_steady_state = 1.24; 
T_steady_state = P_steady_state * sum(R); 

T = 25*ones(n, timesteps);
T(:, 1) = T_steady_state;

% Crank-Nicolson method
for t = 2:timesteps-1
    for j = 1:n-1
        if j == 1
           f = P_steady_state + T(j+1, t) - T(j, t) / (C(j) * R(j));
        else
           f = (T(j-1, t) - T(j, t))/(C(j) * R(j-1)) + (T(j+1, t) - T(j, t))/(C(j) * R(j));
        end
        if j == 1
           f_t = P_steady_state + T(j+1, t+1) - T(j, t+1) / (C(j) * R(j));
        else
           f_t = (T(j-1, t+1) - T(j, t+1))/(C(j) * R(j-1)) + (T(j+1, t+1) - T(j, t+1))/(C(j) * R(j));
        end
        T(j, t+1) = T(j, t) + 0.5 * dt * f_t + 0.5 * dt * f;
    end
end

[Manikanta Aditya]

Hi,

Check the code that implements the Crank-Nicolson method and the Newton-Raphson iterative solver for the transient temperature problem you described:

% Problem parameters
n = 4; % Number of nodes
dt = 0.01; % Time step size
timesteps = 1000; % Number of time steps
R = [6.9364, 8.3004, 2.8952, 12.0162]; % Thermal resistances
C = [1, 2, 1, 2]; % Thermal capacities
P_steady_state = 1.24; % Steady-state power input

% Initial conditions
T_steady_state = P_steady_state * sum(R); % Steady-state temperature
T = T_steady_state * ones(n, 1); % Initialize temperature vector

% Time loop
for t = 2:timesteps
    
    % Newton-Raphson iteration
    max_iter = 100; % Maximum number of iterations
    tol = 1e-6; % Convergence tolerance
    for iter = 1:max_iter
        
        % Compute residuals
        F = zeros(n, 1);
        F(1) = P_steady_state + (T(2) - T(1)) / (C(1) * R(1)) - (dt / 2) * F(1);
        for i = 2:n-1
            F(i) = (T(i-1) - T(i)) / (C(i) * R(i-1)) + (T(i+1) - T(i)) / (C(i) * R(i)) - (dt / 2) * F(i);
        end
        F(n) = (T(n-1) - T(n)) / (C(n) * R(n-1)) - (dt / 2) * F(n);
        
        % Compute Jacobian
        J = zeros(n, n);
        J(1, 1) = -1 / (C(1) * R(1)) - dt / 2;
        J(1, 2) = 1 / (C(1) * R(1));
        for i = 2:n-1
            J(i, i-1) = 1 / (C(i) * R(i-1));
            J(i, i) = -1 / (C(i) * R(i-1)) - 1 / (C(i) * R(i)) - dt / 2;
            J(i, i+1) = 1 / (C(i) * R(i));
        end
        J(n, n-1) = 1 / (C(n) * R(n-1));
        J(n, n) = -1 / (C(n) * R(n-1)) - dt / 2;
        
        % Solve for temperature update
        dT = -J \ F;
        
        % Update temperature
        T = T + dT;
        
        % Check convergence
        if norm(F, Inf) < tol
            break;
        end
    end
    
    % Store temperature for current time step
    T_sol(:, t) = T;
end

Thanks

[larry liu]

@Manikanta Aditya Thank you very much for your help. I change a little code here

% Initial conditions
T_steady_state = 25*ones(n);
for i = 1:n
    for k = i:n
        T_steady_state(i) = T_steady_state(i) + P_steady_state * R(k);
    end
end
T = T_steady_state(:,1) .* ones(n, 1)

If I have a boundary condition where T of each node will be 25 degrees C after steady state, how should I modify it?

The temperature profile of node will be like something like this: