In: Advanced Math
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^'=lny,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
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 %%%%%%%%%%%%%%%%%