Runge kutta giving complex numbers (2024)

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;

Ran in: