The problem is to calculate unknowns and plot graph for Nu vs time.

Matlab code is giving error on odearguments and ode45.

% Main script
xi = 0.1; 
Ra = 100;
R = Ra * xi; 
epsilon = 0.1;

initial_A1_1_1 = 0.01;
initial_A1_1_2 = 0.01;
initial_A1_1_3 = 0.01; 
initial_A1_1_4 = 0.01; 

A0 = [initial_A1_1_1; initial_A1_1_2; initial_A1_1_3; initial_A1_1_4];

% Time span for the simulation
tspan = [0 10]; 

% Solve the ODEs using ode45
[t, A] = ode45(@(t, A) odesystem(t, A, epsilon, Ra, R), tspan, A0);

% Extract the A1_1_1 values from the solution
A1_1_1_values = A(:, 1); 

% Calculate Nu values based on the A1_1_1 values from the ODE solution
Nu_values = xi - 1/2 - (pi^3 / Ra) * A1_1_1_values;

% Plot the results
figure;
plot(t, Nu_values, 'LineWidth', 2);
xlabel('Time');
ylabel('Nu Values');
title('Nu Values vs. Time');
grid on;

function dA = odesystem(t, A, epsilon, Ra, R)
    A1_1_1 = A(1);
    A1_1_2 = A(2);
    A1_1_3 = A(3);
    A1_1_4 = A(4);
    
    % ODEs based on the given expressions
    dA1_1_1 = 3 * epsilon * A1_1_1 - A1_1_3;
    dA1_1_2 = A1_1_2 * A1_1_3 - 100;
    dA1_1_3 = A1_1_3 + Ra * A1_1_4 * A1_1_1;
    dA1_1_4 = A1_1_1;
    
    % Return derivatives as a column vector
    dA = [dA1_1_1; dA1_1_2; dA1_1_3; dA1_1_4];
end

% Call the main function to execute the script

main

Answer 1

[nick]

Hi Madhvi,

I understand that you want to compute the unknowns and plot the graph of Nu vs time.

It seems in the code shared above that the 'main' function is not defined which leads to an error. Additionally, I noticed that some of the 'Nu_values' are resulting in 'NaN' or 'Inf'. This issue arises because the output 'A' from the 'ode45' function contains 'NaN' values or values that are excessively large.

To resolve this, you have to verify the function handle 'odesystem' that defines the system of differential equations to be integrated, as well as the time interval 'tspan' used for the integration are finite.

You can refer to the following documentation to learn more about 'ode45':

https://www.mathworks.com/help/matlab/ref/ode45.html

If you have any further questions, please share the desired output and the derivative of the function to be integrated along with the query.

Hope this helps.

Answer 2

[Sam Chak]

Hi @MADHVI

Your original code runs without error in MATLAB online. However, I have included annotations in the comments to explain why the state variables are unstable. I attempted to address the stability issue to ensure that the Nu variable becomes stable. However, the state variable A1_1_2 remains unstable; nonetheless, this does not affect Nu.

Please note that the Nu variable becomes stable only through my intervention that "metaphysically" alters the original dynamics. These altered dynamics may lack a basis in reality. So, I recommend that you review your differential equations clearly.

% Solve the ODEs using ode45
tspan   = [0 10]; 
initial_A1_1_1 = 0.01;
initial_A1_1_2 = 0.01;
initial_A1_1_3 = 0.01; 
initial_A1_1_4 = 0.01; 
A0      = [initial_A1_1_1; initial_A1_1_2; initial_A1_1_3; initial_A1_1_4];
[t, A]  = ode45(@samchakfix, tspan, A0);

% Calculate Nu values based on the A1_1_1 values from the ODE solution
xi      = 0.1; 
Ra      = 100;
Nu_values = xi - 1/2 - (pi^3 / Ra)*A(:, 1);

% Plot the results
figure
plot(t, A(:,[1, 3, 4]), 'LineWidth', 2), grid on, xlabel('Time / s')
legend('A_1', 'A_3', 'A_4')

figure;
plot(t, Nu_values, 'LineWidth', 2), grid on
xlabel('Time');
ylabel('Nu Values');
title('Nu Values vs. Time');
grid on;

function dA = samchakfix(t, A)
    % parameters
    xi      = 0.1; 
    Ra      = 100;
    R       = Ra * xi; 
    epsilon = 0.1;
    
    % definitions
    A1_1_1  = A(1);
    A1_1_2  = A(2);
    A1_1_3  = A(3);
    A1_1_4  = A(4);
    
    % ODEs based on the given expressions
    % dA1_1_1 = 3*epsilon*A1_1_1 - A1_1_3;    % unstable due to A1_1_1 term
    % dA1_1_2 = A1_1_2*A1_1_3 - 100;          % unstable due to A1_1_3 term
    % dA1_1_3 = A1_1_3 + Ra*A1_1_4*A1_1_1;    % unstable due to A1_1_3 term
    % dA1_1_4 = A1_1_1;                       % unstable since A1_1_1 is unstable
    
    % My fix
    dA1_1_1 = - 3*epsilon*A1_1_1 - A1_1_3;  % stability depends on A1_1_3 term
    dA1_1_2 = A1_1_2*1*A1_1_3 - 100;        % unstable due to A1_1_3 term, but doesn't affect Nu
    dA1_1_3 = - A1_1_3 + Ra*A1_1_4*A1_1_1;  % stability depends on A1_1_4 term
    dA1_1_4 = A1_1_1;                       % stability depends on A1_1_1 term
    
    % Return derivatives as a column vector
    dA      = [dA1_1_1; 
               dA1_1_2; 
               dA1_1_3; 
               dA1_1_4];
end