Index in position 1 is invalid. Array indices must be positive integers or logical values.

5 views (last 30 days)
function mutind2()
h = 0.0001;
t = 0;
T = [];
y_euler = [0; 0; 0];
y_heun = [0; 0; 0];
y_rk4 = [0; 0; 0];
e_euler = [];
e_heun = [];
e_rk4 = [];
f = 5; % Frequency in Hz
while t < 30
E = 210 * sin(2*pi*f*t);
y_euler = [y_euler, y_euler(:, end) + h * f(t, y_euler(:, end), E(:, end))];
y_heun = [y_heun, heun_step(t, y_heun(:, end), h, E(:, end))];
y_rk4 = [y_rk4, rk4_step(t, y_rk4(:, end), h, E(:, end))];
e_euler(end+1) = E;
e_heun(end+1) = E;
e_rk4(end+1) = E;
t = t + h;
T(end+1) = t;
end
% Plotting currents (i1 and i2) for each excitation with Euler method
figure(1);
plot(T, y_euler(1, 2:end), 'r', T, y_euler(2, 2:end), 'b');
legend('i1 Euler', 'i2 Euler');
xlabel('t [s]');
ylabel('i [A]');
grid on;
title('Currents (Euler Method)');
% Plotting currents (i1 and i2) for each excitation with Heun's method
figure(2);
plot(T, y_heun(1, 2:end), 'r', T, y_heun(2, 2:end), 'b');
legend('i1 Heun', 'i2 Heun');
xlabel('t [s]');
ylabel('i [A]');
grid on;
title("Currents (Heun's Method)");
% Plotting currents (i1 and i2) for each excitation with RK4 method
figure(3);
plot(T, y_rk4(1, 2:end), 'r', T, y_rk4(2, 2:end), 'b');
legend('i1 RK4', 'i2 RK4');
xlabel('t [s]');
ylabel('i [A]');
grid on;
title('Currents (RK4 Method)');
% Plotting voltages (uc and e) for each excitation with Euler method
figure(4);
plot(T, y_euler(3, 2:end), 'r', T, e_euler, 'b');
legend('uc Euler', 'e Euler');
xlabel('t [s]');
ylabel('Voltage');
grid on;
title('Voltages (Euler Method)');
% Plotting voltages (uc and e) for each excitation with Heun's method
figure(5);
plot(T, y_heun(3, 2:end), 'r', T, e_heun, 'b');
legend('uc Heun', 'e Heun');
xlabel('t [s]');
ylabel('Voltage');
grid on;
title("Voltages (Heun's Method)");
% Plotting voltages (uc and e) for each excitation with RK4 method
figure(6);
plot(T, y_rk4(3, 2:end), 'r', T, e_rk4, 'b');
legend('uc RK4', 'e RK4');
xlabel('t [s]');
ylabel('Voltage');
grid on;
title('Voltages (RK4 Method)');
disp('Simulation results:');
disp('-------------------');
disp(['Final value of U-C: ' num2str(y_euler(3, end)) ' V']);
disp(['Final value of i-1: ' num2str(y_euler(1, end)) ' A']);
disp(['Final value of i-2: ' num2str(y_euler(2, end)) ' A']);
end
function dy = f(t, y, E)
R1 = 0.1;
R2 = 10;
C = 0.5;
L1 = 3;
L2 = 5;
M = 0.8;
D1 = L1 / M - M / L2;
D2 = M / L1 - L2 / M;
A = [-R1 / (M * D1), R2 / (L2 * D1), -1 / (M * D1);
-R1 / (L1 * D2), R2 / (M * D2), -1 / (L1 * D2);
1 / C, 0, 0];
b = [1 / (M * D1); 1 / (L1 * D2); 0];
dy = A * y + b * E;
end
function y_next = euler_step(t, y, h, E)
y_next = y + h * f(t, y, E);
end
function y_next = heun_step(t, y, h, E)
k1 = f(t, y, E);
k2 = f(t + h, y + h * k1, E);
y_next = y + (h / 2) * (k1 + k2);
end
function y_next = rk4_step(t, y, h, E)
k1 = f(t, y, E);
k2 = f(t + h/2, y + (h/2) * k1, E);
k3 = f(t + h/2, y + (h/2) * k2, E);
k4 = f(t + h, y + h * k3, E);
y_next = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4);
end

Answers (2)

Sahas Marwah
Sahas Marwah on 12 Jun 2023
Hi @Dmytro,
After running your code, I found an error in line 14.
You have defined f = 5 (in line 11) and used it as a function (in line 14). Hence, it gives out the error.
I am unable to understand how you want to use f = 5 at line 14. So, please correct this line according to your functionality.

Image Analyst
Image Analyst on 12 Jun 2023
It doesn't know whether, when you go to assign y, whether you want to call your function f or want to use the f you assigned a few lines above.
It's thoroughly explained in the FAQ:
Another thing that's wrong is there are far too few comments. For example the functions don't even describe what they do.

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!