Runge kutta giving complex numbers – urgent please :)
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:TESIDati29gen_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;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:TESIDati29gen_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; 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:TESIDati29gen_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; differential equations, complex numbers MATLAB Answers — New Questions
