Runge kutta giving complex numbers

Question

SIMONE GUAZZOTTI 2024-6-27，18:16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2132656-runge-kutta-giving-complex-numbers

评论： Umar 2024-6-27，23:10

29gen_1feb_2023_z_hub.csv

在 MATLAB Online 中打开

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;

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Umar 2024-6-27，18:31

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2132656-runge-kutta-giving-complex-numbers#answer_1478066

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

2 个评论
显示无隐藏无

SIMONE GUAZZOTTI 2024-6-27，22:18

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 2024-6-27，23:10

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

请先登录，再进行评论。

Answer 2

Torsten 2024-6-27，20:01

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2132656-runge-kutta-giving-complex-numbers#answer_1478096

编辑：Torsten 2024-6-27，21:30

在 MATLAB Online 中打开

29gen_1feb_2023_z_hub.csv

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Runge kutta giving complex numbers

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

2 个评论
显示无隐藏无

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Runge kutta giving complex numbers

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

2 个评论 显示 无隐藏 无

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

WeChat

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

2 个评论
显示无隐藏无

0 个评论
显示 -2更早的评论隐藏 -2更早的评论