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

13 views (last 30 days)
Dmytro on 12 Jun 2023
Answered: Image Analyst on 12 Jun 2023
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
##### 0 CommentsShow -2 older commentsHide -2 older comments

Sign in to comment.

### Answers (2)

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.
##### 0 CommentsShow -2 older commentsHide -2 older comments

Sign in to comment.

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.
##### 0 CommentsShow -2 older commentsHide -2 older comments

Sign in to comment.

### 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!