Question

The algorithm is basically as follows. The notation is slightly different from that in the website you were given, but there is no difference in the method.

Given the initial value problem

dy/dx=f(x,y),y(a)= y_0

Euler’s Method with step size h consists in applying the iterative formula

y_(n+1)= y_n+h∙f(x_n,y_n ),n≥0

To compute successive approximations y_1,y_2,y_3,⋯ to the (true) values 〖y(x〗_1),〖y(x〗_2),〖y(x〗_3),⋯ of the exact solution y=y(x) at the points x_1,x_2,x_3,⋯, respectively.

In plain English:
You want to approximate the value of dy/dx (or y’) at some point in an interval.

Step 1: Depending on how accurate you need to be, divide the interval up into little pieces of equal length; this length is the step size h. For purposes of discussion, let’s use the interval [0,1] and use ten intervals, so h = 0.1.

Step 2: y_0=0
Step 3: y_1=y_0+0.1f(x_0,y_0)
Step 4: y_2=y_1+0.1f(x_1,y_1)
…
Stop after ten steps, in this case. Usually the stopping criterion is a level of accuracy.

You can easily set this up in Excel.

Exercises
Use Euler’s Method with step sizes h =0.1,0.02, 0.004, 0.0008 (that is, do the problem 4 times, each with a more precise value of h) , 10 equally spaced iterations.


1. y^'=x^2+y^2,y(0)=0,0≤x≤1

2. y^'=x^2-y^2,y(0)=1,0≤x≤2

3. y^'=ln⁡y,y(1)=2,1≤x≤2

4. y^'=x^(2/3)+y^(2/3),y(0)=1,0≤x≤2

5. y^'=x+√x,y(0)=1,0≤x≤2

6. y^'=x+∛x,y(0)= -1,0≤x≤2

Answer 1

clear all
close all

%function for which Euler method have to calculate
f=@(y,x) x.^2+y.^2;
%all h
hh=[0.1 0.02 0.004 0.0008];
fprintf('For the function f(x,y)=\n')
disp(f)

for i=1:length(hh)
    %step size h
    h=hh(i);
    %initial value
    y0=0; x0=0;
    x_end=1;
    [y1_result,x_result] = euler_method(f,y0,x0,x_end,h);
    figure(1)
    hold on
    plot(x_result,y1_result)
    lgnd{i}=sprintf('For step size %.4f',h);
end
xlabel('x')
ylabel('y(x)')
title('y(x) vs. x plot for function 1')
legend(lgnd)
box on

f=@(y,x) x.^2-y.^2;
%all h
hh=[0.1 0.02 0.004 0.0008];
fprintf('For the function f(x,y)=\n')
disp(f)

for i=1:length(hh)
    %step size h
    h=hh(i);
    %initial value
    y0=1; x0=0;
    x_end=2;
    [y1_result,x_result] = euler_method(f,y0,x0,x_end,h);
    figure(2)
    hold on
    plot(x_result,y1_result)
    lgnd{i}=sprintf('For step size %.4f',h);
end
xlabel('x')
ylabel('y(x)')
title('y(x) vs. x plot for function 2')
legend(lgnd)
box on

f=@(y,x) log(y);
%all h
hh=[0.1 0.02 0.004 0.0008];
fprintf('For the function f(x,y)=\n')
disp(f)

for i=1:length(hh)
    %step size h
    h=hh(i);
    %initial value
    y0=2; x0=1;
    x_end=2;
    [y1_result,x_result] = euler_method(f,y0,x0,x_end,h);
    figure(3)
    hold on
    plot(x_result,y1_result)
    lgnd{i}=sprintf('For step size %.4f',h);
end
xlabel('x')
ylabel('y(x)')
title('y(x) vs. x plot for function 3')
legend(lgnd)
box on

f=@(y,x) x.^(2/3)+y.^(2/3);
%all h
hh=[0.1 0.02 0.004 0.0008];
fprintf('For the function f(x,y)=\n')
disp(f)

for i=1:length(hh)
    %step size h
    h=hh(i);
    %initial value
    y0=1; x0=0;
    x_end=2;
    [y1_result,x_result] = euler_method(f,y0,x0,x_end,h);
    figure(4)
    hold on
    plot(x_result,y1_result)
    lgnd{i}=sprintf('For step size %.4f',h);
end
xlabel('x')
ylabel('y(x)')
title('y(x) vs. x plot for function 4')
legend(lgnd)
box on

f=@(y,x) x+sqrt(x);
%all h
hh=[0.1 0.02 0.004 0.0008];
fprintf('For the function f(x,y)=\n')
disp(f)

for i=1:length(hh)
    %step size h
    h=hh(i);
    %initial value
    y0=1; x0=0;
    x_end=2;
    [y1_result,x_result] = euler_method(f,y0,x0,x_end,h);
    figure(5)
    hold on
    plot(x_result,y1_result)
    lgnd{i}=sprintf('For step size %.4f',h);
end
xlabel('x')
ylabel('y(x)')
title('y(x) vs. x plot for function 5')
legend(lgnd)
box on

f=@(y,x) x+x.^(1/3);
%all h
hh=[0.1 0.02 0.004 0.0008];
fprintf('For the function f(x,y)=\n')
disp(f)

for i=1:length(hh)
    %step size h
    h=hh(i);
    %initial value
    y0=-1; x0=0;
    x_end=2;
    [y1_result,x_result] = euler_method(f,y0,x0,x_end,h);
    figure(6)
    hold on
    plot(x_result,y1_result)
    lgnd{i}=sprintf('For step size %.4f',h);
end
xlabel('x')
ylabel('y(x)')
title('y(x) vs. x plot for function 6')
legend(lgnd)
box on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Matlab code for Euler's forward
function [y1_result,t_result] = euler_method(f,y0,t0,tend,h)
%function for Euler equation solution
  
    %all step size
    N=(tend-t0)/h;
    %Initial values

    %t end values
    tn=t0:h:tend;
    % Euler steps
    y1_result(1)=y0;
    t_result(1)=t0;
  
    for i=1:length(tn)-1
        t_result(i+1)= t_result(i)+h;
        y1_result(i+1)=y1_result(i)+h*double(f(y1_result(i),t_result(i)));
    end
end

%%%%%%%%%%%%%%%%%%% End of Code %%%%%%%%%%%%%%%%%