Forward Euler solution improvement

I would like to know if this code can be improved or is well-written as it is? it does not give errors. I'll appreciate any comments.Thanks

%part a

clear; clc;

% Parameters
alpha = 10;
beta_values = [2, 4];
h = 0.01;
t = 0:h:20;
N = length(t);

% Initial conditions
y1_0 = 0;
y2_0 = 2;

for b = 1:length(beta_values)
    beta = beta_values(b);

    
    y1 = zeros(1, N);
    y2 = zeros(1, N);

    % Initial values
    y1(1) = y1_0;
    y2(1) = y2_0;

    % Forward Euler method
    for n = 1:N-1
        f1 = alpha ...
             - y1(n) ...
             - (4*y1(n)*y2(n)) / (1 + y1(n)^2);

        f2 = beta * y1(n) * ...
             (1 - y2(n) / (1 + y1(n)^2));

        y1(n+1) = y1(n) + h * f1;
        y2(n+1) = y2(n) + h * f2;
    end

    figure;
    plot(t, y1, 'LineWidth', 1.5);
    xlabel('t');
    ylabel('y_1(t)');
    title(['Forward Euler: y_1 vs t, \beta = ', num2str(beta)]);
    grid on;

    figure;
    plot(y1, y2, 'LineWidth', 1.5);
    xlabel('y_1');
    ylabel('y_2');
    title(['Phase Portrait: y_2 vs y_1, \beta = ', num2str(beta)]);
    grid on;
end

%part b


% Parameters
alpha = 10;
beta_values = [3.4, 3.5, 3.6];   
h = 0.01;
t = 0:h:100;                    
N = length(t);

% Initial conditions
y1_0 = 0;
y2_0 = 2;

for b = 1:length(beta_values)
    beta = beta_values(b);

    
    y1 = zeros(1, N);
    y2 = zeros(1, N);

    % Initial values
    y1(1) = y1_0;
    y2(1) = y2_0;

    % Forward Euler method
    for n = 1:N-1
        f1 = alpha ...
             - y1(n) ...
             - (4*y1(n)*y2(n)) / (1 + y1(n)^2);

        f2 = beta * y1(n) * ...
             (1 - y2(n) / (1 + y1(n)^2));

        y1(n+1) = y1(n) + h * f1;
        y2(n+1) = y2(n) + h * f2;
    end

    
    figure;
    plot(t, y1, 'LineWidth', 1.5);
    xlabel('t');
    ylabel('y_1(t)');
    title(['y_1 vs t, \beta = ', num2str(beta)]);
    grid on;

    figure;
    plot(y1, y2, 'LineWidth', 1.5);
    xlabel('y_1');
    ylabel('y_2');
    title(['Phase Portrait y_2 vs y_1, \beta = ', num2str(beta)]);
    grid on;
end

Answer 1

[Alan Stevens]

It's clear and runs quickly. However, it's probably neater to pre-define your gradient functions like this:

% Gradient functions
f1 = @(alpha,y1,y2) alpha - y1 - (4*y1*y2) / (1 + y1^2);
f2 = @(beta,y1,y2)  beta * y1 * (1 - y2 / (1 + y1^2));

and then replace the forward Euler equations simply by:

    % Forward Euler method
    for n = 1:N-1
        y1(n+1) = y1(n) + h * f1(alpha,y1(n),y2(n));
        y2(n+1) = y2(n) + h * f2(beta,y1(n),y2(n));
    end

[Cesar]

Thanks, I would like to know if this code can also be improved? it works and does not give errors. Any comments will be appreciated.

clc; clear; close all;

% Parameters

T_final = 10;

u0 = [1; 0];

N_steps = [10000, 30000, 50000];

figure;

hold on;

colors = ['r', 'g', 'b'];

for i = 1:length(N_steps)

N = N_steps(i);

h = T_final / N;

t = 0:h:T_final;

u = zeros(2, length(t));

u(:, 1) = u0;

% Forward Euler

for n = 1:N

f = [u(2, n); -u(1, n)];

u(:, n+1) = u(:, n) + h * f;

end

plot(u(1, :), u(2, :), 'Color', colors(i), 'LineWidth', 1, ...

'DisplayName', ['N = ' num2str(N)]);

end

xlabel('u_1');

ylabel('u_2');

title('u_2 vs u_1 (Forward Euler)');

set(gcf, 'Position', [100, 100, 00, 500])

legend('show');

grid on;

axis equal;

[Cesar]

Also just wondering if oed45 is correctly implemented here?:

function orbit

mu = 0.012277471;

mu_hat = 1 - mu;

T_final = 17.06521656; % One full period

y0 = [0.994; 0; 0; -2.00158510637908252240537862224];

options = odeset('RelTol', 1e-12, 'AbsTol', 1e-12);

[t, y] = ode45(@(t, y) arenstorf_system(t, y, mu, mu_hat), [0 T_final], y0, options);

figure('Color', 'w', 'Position', [100, 100, 900, 700]);

plot(y(:,1), y(:,2), 'b--', 'LineWidth', 1.5);

hold on;

xlabel('u_1'); ylabel('u_2');

legend('ode45 solution');

grid on;

axis equal; % Ensures the orbit isn't visually distorted

end

function dy = arenstorf_system(~, y, mu, mu_hat)

u1 = y(1); u2 = y(2); v1 = y(3); v2 = y(4);

D1 = ((u1 + mu)^2 + u2^2)^(3/2);

D2 = ((u1 - mu_hat)^2 + u2^2)^(3/2);

du1 = v1;

du2 = v2;

dv1 = u1 + 2*v2 - mu_hat*(u1 + mu)/D1 - mu*(u1 - mu_hat)/D2;

dv2 = u2 - 2*v1 - mu_hat*u2/D1 - mu*u2/D2;

dy = [du1; du2; dv1; dv2];

end

[Alan Stevens]

Your code posted at 19:56 above only puts titles, labels etc on the last plot (as the circles are all the same size, you only see the blue one). Also the set gcf command confuses the axis equal command. You could do the following:

T_final = 10; 
u0 = [1; 0];  
N_steps = [10000, 30000, 50000];

colors = ['r', 'g', 'b'];
for i = 1:length(N_steps)
    N = N_steps(i);
    h = T_final / N;
    t = 0:h:T_final;
    
    u = zeros(2, length(t));
    u(:, 1) = u0;
    
    % Forward Euler
    for n = 1:N
       
        f = [u(2, n); -u(1, n)]; % du1/dt = u2   du2/dt = -u1
        
        u(:, n+1) = u(:, n) + h * f;
    end
    %figure  - if you want three separate figures
    subplot(2,2,i) % - if you want all tyhree on one figure
    plot(u(1, :), u(2, :), 'Color', colors(i), 'LineWidth', 1, ...
        'DisplayName', ['N = ' num2str(N)]);
    xlabel('u_1');
    ylabel('u_2');
    title('u_2 vs u_1 (Forward Euler)');
    %set(gcf, 'Position', [100, 100, 00, 500])
    legend('show');
    grid on;
    axis equal;
end

[Alan Stevens]

Your ode45 is correctly implemented. (The fourth value of y0 seems to be unneccesarily precise! However, I didn't test the sensitivity of the results to reducing the number of decimal places.)