Forward Euler solution improvement

Question

Cesar 2026-2-6，1:46

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2182524-forward-euler-solution-improvement

评论： Cesar about 8 hours 前

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

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

请先登录，再进行评论。

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

Answer 1

Alan Stevens about 19 hours 前

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2182524-forward-euler-solution-improvement#answer_1573645

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

Cesar about 2 hours 前

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 about 8 hours 前

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

请先登录，再进行评论。

Forward Euler solution improvement

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

采纳的回答

2 个评论
显示无隐藏无

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Forward Euler solution improvement

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

采纳的回答

2 个评论 显示 无隐藏 无

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

2 个评论
显示无隐藏无