Runge kutta giving complex numbers

Hello everyone,

I need to graduate in really few days. Please if you could help me I would really appreciate it.

I need to solve a differential equation in matlab, but during the simulation V1_rel becomes a complex number ( in particular at the timestep t1(7617) where V1_rel(7617) = V1(7617) + v1_nac(7616) , but these last two values are real values). So somehow a sum of two real values result in a complex number. From that iteration onwards also v1_nac become complex (so from v1_nac(7617).

But the problem probably starts before in the calculation of v1, because it starts low but rapidly increases too much in absolute value. v1 is in rad/s and it must be below 0.01 . Consequently v1_nac is too high and V1_rel assumes also negative values, which are not phisically possible.

I add at the bottom the script with which I derived the V1 values, but surely they are not the problem. They have been used for other simulations without problems. I attached also the file from which I derived the V1 values.

function dydt = odefun(t, yv, I, B, K, M)
    y = yv(1);
    v = yv(2);
    dy_dt = v;
    dv_dt = M/I - (B/I) * v  - K/I * y - M/I * y^2;
    dydt = [dy_dt; dv_dt];
end

%% Δt = 1s
I = 81691700806.8323;
B = 247901197.414257;
K = 1880699112.40857;

y1 = zeros(length(t1), 1);               
y1_deg = zeros(length(t1),1);             
v1 = zeros(length(t1), 1);               
v1_nac = zeros(length(t1), 1);            
V1_rel = zeros(length(t1), 1);           
F1_thrust = zeros(length(t1),1);         
M1_thrust = zeros(length(t1),1);

y1(1) = 0.005;
v1(1) = 0.001;
v1_nac(1) = 0.01;
V1_rel(1) = V1(1);
F1_thrust(1) = (6.073 * V1(1)^1.9808) *1000;
M1_thrust(1) = F1_thrust(1) * 105;
% Metodo di Runge-Kutta di quarto ordine
for i = 2:length(t1)
    dt = t1(i) - t1(i-1);
   
    V1_rel(i)= V1(i) + v1_nac(i-1);
    
    if  V1_rel(i) < 11
        F1_thrust(i) =  (6.073 *  power( V1_rel(i),1.9808) ) *1000 ;
    else
        F1_thrust(i) = (-0.0455 *  power(V1_rel(i),4) + 3.2312* power(V1_rel(i),3) - 84.841* power(V1_rel(i),2) + 943.93* V1_rel(i) - 3036.5) * 1000  ;
    end 
    M1_thrust(i) = F1_thrust(i) * 105;
    M = M1_thrust(i);
    yv1 = [y1(i-1); v1(i-1)];
   
    k1_ = odefun(t1(i-1), yv1, I, B, K, M);
    k2_ = odefun(t1(i-1) + dt/2, yv1 + dt/2 * k1_, I, B, K, M);
    k3_ = odefun(t1(i-1) + dt/2, yv1 + dt/2 * k2_, I, B, K, M);
    k4_ = odefun(t1(i), yv1 + dt * k3_, I, B, K, M);
    
    % Aggiornamento delle soluzioni
    yv_new = yv1 + dt/6 * (k1_ + 2*k2_ + 2*k3_ + k4_);
    y1(i) =  yv_new(1) ;  
    v1(i) =  yv_new(2) ;  
    v1_nac(i) = v1(i)*105;    % from rad/10s to m/s%
    y1_deg(i) = rad2deg(y1(i));
end

% Read the data from the CSV file. Extract time and velocity
data = readtable('C:\TESI\Dati\29gen_1feb_2023_z_hub.csv');
time = data{:, 1};
velocity = data{:, 2};

% Convert time to datetime format
time = datetime(time, 'InputFormat', 'dd/MM/yyyy HH:mm:ss');
time_seconds = seconds(time - time(1));

% Interpolate linearly velocity every 10 seconds
t10 = (0:10:max(time_seconds))';
t1= (time_seconds(338):1:time_seconds(386))';  % 30 jan 09:00 - 30 jan 18:00
V_interp10 = interp1(time_seconds, velocity, t10, 'linear');
V_interp1 = interp1(time_seconds, velocity, t1, 'linear');

% Add noise proportional to wind variability if needed
n = randn(length(t10), 1) / 35;
vett2 = [diff(V_interp10); 0];
%vett2_norm = ( vett2 / std(vett2) );
%vett2_adjusted = 0.95 * vett2_norm + (1 - 0.95) * 1;
%n = n2 .* vett2_adjusted * 7
n10 = n .* ( vett2 / std(vett2) ) * 7;       
n1 = interp1(t10, n10 , t1);
                                                                             
V10 = V_interp10 + n10;                       
V1 = V_interp1 + n1;

Answer 1

[Umar]

Hi Simone,

After reviewing your code, the problem lies in the calculation of v1. It seems that the formula used to update v1 is not correctly capturing the desired behavior. As a result, v1 increases too rapidly, leading to the subsequent issues with v1_nac and V1_rel. To fix the issue, modify the formula used to update v1 in the code. Currently, the formula is:

dv_dt = M/I - (B/I) * v - K/I * y - M/I * y^2;

This formula is causing v1 to increase too rapidly. We can try adjusting the coefficients to slow down the increase of v1. For example, we can divide the coefficient (B/I) by a larger value, such as 10, to reduce its impact on the increase of v1. The modified formula would be:

dv_dt = M/I - (B/I)/10 * v - K/I * y - M/I * y^2;

By making this adjustment, we can control the increase of v1 and prevent it from reaching excessively high values.

It sounds like you are trying to model the thrust and velocity of a system using Matlab and using the Runge-Kutta method to solve the ODEs that describe the system's behavior. Interesting

[SIMONE GUAZZOTTI]

Hello Umar,

Thank you very much for your answer. Indeed, I should adjust the (B/I) term and dividing it for 0.1 it gives reasonable and possible values. I guess the only I have to find the optimal values is to try different values and choose according the result.

Anyways by the way, y and v represent the inclination angle and velocity of oscillation of a floating wind turbine.

Thank you very much for your availability and suggestion.

[Umar]

Glad to help out, Simone. Kudos to Torsten who contributed efforts to resolve this issue as well.

Answer 2

[Torsten]

It doesn't make sense to use such discontinuous input data for the solution of a differential equation.

Solving a differential equation (numerically) means: the solution function is differentiable. For this, the inputs should be at least continuous. Not the case here.

% Read the data from the CSV file. Extract time and velocity
data = readtable('29gen_1feb_2023_z_hub.csv');
time = data{:, 1};
velocity = data{:, 2};
% Convert time to datetime format
time = datetime(time, 'InputFormat', 'dd/MM/yyyy HH:mm:ss');
time_seconds = seconds(time - time(1));
hold on
plot(time_seconds,velocity)
% Interpolate linearly velocity every 10 seconds
t10 = (0:10:max(time_seconds))';
t1= (time_seconds(338):1:time_seconds(386))';  % 30 jan 09:00 - 30 jan 18:00
V_interp10 = interp1(time_seconds, velocity, t10, 'linear');
V_interp1 = interp1(time_seconds, velocity, t1, 'linear');
% Add noise proportional to wind variability if needed
n = randn(length(t10), 1) / 35;
vett2 = [diff(V_interp10); 0];
%vett2_norm = ( vett2 / std(vett2) );
%vett2_adjusted = 0.95 * vett2_norm + (1 - 0.95) * 1;
%n = n2 .* vett2_adjusted * 7
n10 = n .* ( vett2 / std(vett2) ) * 7;       
n1 = interp1(t10, n10 , t1);
                                                                             
V10 = V_interp10 + n10;                       
V1 = V_interp1 + n1;
plot(t1,V1)
hold off