I need to solve a cauchy problem with the Euler Method.

Hi all. I was asked to solve this problem by my teacher:

I have to write a function that solves this cauchy problem with the Eulerian method, using an h (step size) of 0.25, in the interval [0,2].

I am then told:

"knowing that the exact solution of the equation is y(x)=cos(x)+sin(x), calculate the error committed in estimating y(2)."

I have to also make a graph of both the true solution (in one colour) and the solution made with the eulerian method in another colour.

What i've done so far (Before reaching the limit of my capabilities!):

I've written this script:

% Eulerian Method
% Setup and initial conditions
h = 0.25;  % step size
x = 0:h:5;  % interval of x
y = zeros(size(x));  % allocate the resulting y
y(1) = 1;  % inizial value of y
n = numel(y);  % the number of y values
% The loop that solves the differential equation:
for i=1:n-1
    f = (y(i)-2*sin(x(i)))
    y(i+1) = y(i) + h * f;
end

And i've managed to plot it into a graph but i have no idea of how to complete the tasks properly.

Please Help!

Answer 1

[Ameer Hamza]

Following shows how to find true solution using symbolic toolbox and plot both of them on same graph

% Eulerian Method
% Setup and initial conditions
h = 0.25;  % step size
x = 0:h:5;  % interval of x
y = zeros(size(x));  % allocate the resulting y
y(1) = 1;  % inizial value of y
n = numel(y);  % the number of y values
% The loop that solves the differential equation:
for i=1:n-1
    f = (y(i)-2*sin(x(i)));
    y(i+1) = y(i) + h * f;
end

syms Y(X)
eq = diff(Y) == Y - 2*sin(X);
ic = Y(0) == 1;
sol = dsolve(eq, ic);
y_sol = double(subs(sol, X, x));

hold on
plot(x, y_sol, 'r', 'DisplayName', 'true solution');
plot(x, y, 'b', 'DisplayName', 'Euler method');
legend

[Stefano Baioli]

Thank you for your help!

At the end of the script, i've added this, in order to calculate the error:

y1=(interp1(x,y,2))
y2=(interp1(x,y_sol,2))
Error= y1-y2 %Calculates the error

Do you think this is all the teacher was asking of me? Thank you!

[Ameer Hamza]

Yes, the error estimation seems correct, according to the question statement.

[Edinson Manga]

Hello! and if I want to have several h how would I change the algorithm?

[John D'Errico]

I would suggest it is a bad idea to use interp1 there, certainly so if you use the default method for interp1. The default is LINEAR interpolation. In that case, on a coarse interval, a significant fraction of the error you might report could just be that of interp1!

In context of what you are doing, it just so happens that x==2 just happens to be one of the points where you predicted the solution. In that case, linear interpolation is meaningless, since you would just want to see the value at the point of interest, and interp1 will do exactly the same. However, suppose the question had been to predict the error at some intermediate point? In that case curvature of the result would be a bad thing, since the connect-the-dots linear interpolation done inside interp1 will hurt you.

In the specific case here where you just need the error at x==2, you could have done just this:

Error = y(9) - (sin(x) + cos(x));

In a more general case, perhaps where you were asked to provide the error at say x=2.13, thus intermediate to your data, now you might more accurately have done:

xerr = 2.13;
ytrue = @(x) sin(x) + cos(x);
Error = interp1(x,y,xerr,'spline') - ytrue(xerr);

At least here the error of interpolation will be reduced due to the use of a spline, and only for the vector y. See here that I used the true function to predict the value at the point of interest, not an interpolation at all.

Always try to avoid adding additional sources of error to what you do.

An even better choice would have been to carefully write way-points into your code, thus making sure the ODE solver predicts a value at exactly some points of interst. This requires more carefully written code of course.

[Ameer Hamza]

Yes, interp1 will introduce its own error, however, in the previous case, I ignored it because the query point xq=2 also lies in vector x. However, in the general case, it can fail. Following code is for multiple values of h and calculates error using the method suggested by John

% Eulerian Method
% Setup and initial conditions
H = [0.25 0.1 0.05 0.01 0.005];  % step size
Y_euler = cell(1,numel(H));
for k=1:numel(H)
    h = H(k);
    x = 0:h:5;  % interval of x
    y = zeros(1, numel(x));  % allocate the resulting y
    y(1) = 1;  % inizial value of y
    n = numel(y);  % the number of y values
    % The loop that solves the differential equation:
    for i=1:n-1
        f = (y(i)-2*sin(x(i)));
        y(i+1) = y(i) + h * f;
    end
    Y_euler{k} = [x; y].';
end

syms Y(X)
eq = diff(Y) == Y - 2*sin(X);
ic = Y(0) == 1;
sol = dsolve(eq, ic);
y_sol = double(subs(sol, X, x));

hold on
plot(x, y_sol, 'r', 'DisplayName', 'true solution');
for k=1:numel(H)
    plot(Y_euler{k}(:,1), Y_euler{k}(:,2), 'b--', 'DisplayName', ['Euler method h=' num2str(H(k), '%.3f')]);
end
legend

err = zeros(1,numel(H));
for k=1:numel(H)
    err(k) = interp1(Y_euler{k}(:,1), Y_euler{k}(:,2), 2) - subs(sol, X, 2);
end