Euler Method Setting h Value

Question

I am trying to implement euler method to solve differential equation y' = x^3 + y^2 between 0 and 2 with initial condition y(0) = 0.5.

Firstly I set h = 0.1 and it's okey. y converges to 643,....

Secondly I set h = 0.01 and y diverges to the infinity.

Can anybody explain why this happen?

My MATLAB code is below:

clc;
clear;
close all;

% dy/dt = f(t, y)
f = @(t, y) t.^3 + y.^2;

a = 0;
ya = 0.5;
b = 2;
h = 0.1;
% h = 0.01;
M = (b-a)/h;
E = euler(f, a, b, ya, M);

function [E] = euler(f, a, b, ya, M)
% INPUT  - f is the function
%        - a and b are the left and right endpoints
%        - ya is the initial condition y(a)
%        - M is the number of steps
% OUTPUT - E = [T Y] where T is the vector of abscissas and
%          Y is the vector of ordinates

h = (b-a) / M;

T = (a:h:b)';
Y = zeros(M+1, 1);

Y(1) = ya;

for j = 1:M
   Y(j+1) = Y(j) + h*f(T(j), Y(j));
end

E = [T Y];
end

Lutz Lehmann · Answer 1 · 2020-12-30T08:23:35.160

This happens because the exact solution has a singularity in that interval and with decreasing steps size you more and more approximate the exact solution with its divergence.

You can solve this Riccati equation more safely using the substitution y=-u'/u to get the second order equation u''(x)+x^3*u(x)=0. Roots of u are singularities of y. The code

f = lambda x,u: np.array([u[1], -x**3*u[0]])
u0 = np.array([1,-0.5])
x_ref = np.linspace(0,2,201)
u_ref = odeint(f,u0, x_ref, tfirst=True, atol=1e-8, rtol=1e-10)

x_euler=np.linspace(0,2,21)
dx = x_euler[1]-x_euler[0]
u_euler = np.array(len(x_euler)*[u0])
for i in range(len(x_euler)-1):
    u_euler[i+1] = u_euler[i] + dx*f(x_euler[i],u_euler[i])

plt.plot(x_ref,u_ref[:,0],'b', lw=3)
plt.plot(x_euler, u_euler[:,0], '-xr', lw=1)

gives a plot of

with a root of u_ref at about x=1.5867, and the Euler plot is not too far away.

Euler Method Setting h Value

1 Answers1