# Euler's Method and exact solution plot

15 views (last 30 days)
scott lamb on 26 Nov 2020
Commented: John D'Errico on 26 Nov 2020
when i run the code below, the plots do not look correct to me. As when i make the time step smaller the eulers should get more accurate, thus closer matching the exact solution? But when i do this they eulers and exact solution seem to get further apart on the plot.
clear;
h=0.05; % time step
t=0:h:4; %time range
y=zeros(size(t)); % set y array all 0's to same size as t
y(1)=2; % set initial condition at time 0 to 2
n=numel(y);
dydt = 4*exp(0.8*t) - 0.5*y;
exact_sol=(4/1.3)*(exp(0.8*t)-exp(-0.5*t))+2*exp(-0.5*t); %This is the exact solution to dy/dt
for i=1 : n-1 %for loop to interate through y values for
y(i+1)= y(i)+ h * dydt(i); % the Euler method
end
plot(t,y) %plot Euler
hold on
plot(t,exact_sol,'red'); % plots the exact solution to this differential equation
grid on

John D'Errico on 26 Nov 2020
Edited: John D'Errico on 26 Nov 2020
Sorry, but your code is wrong (so is your exact solution.) Why? Look at what you wrote.
h=0.05; % time step
t=0:h:4; %time range
y=zeros(size(t)); % set y array all 0's to same size as t
y(1)=2; % set initial condition at time 0 to 2
n=numel(y);
dydt = 4*exp(0.8*t) - 0.5*y;
So dydt is a VECTOR. t and y are vectors. However, y is a vector with first element 2, and EVERY other element 0.
for i=1 : n-1 %for loop to interate through y values for
y(i+1)= y(i)+ h * dydt(i); % the Euler method
end
How do you use dydt? Remember, it was fixed before the loop. Those vector elements don't change. You never update dydt when y changes. Therefore you are not solving the problem you think you are solving. You need to update dydt as a function, that changes as both y and t vary. If you don't, well then you get arbitrary garbage as a result - what you got, in fact. That is, you CANNOT create dydt in advance as avector.
So first, let me check your claim as to the solution.
syms y(t)
true_sol = dsolve(diff(y,t,1) == 4*exp(0.8*t) - 0.5*y(t),y(0) == 2)
true_sol = Interestingly, that is also not the same as what you claim to be the solution. So exact_sol is also incorrect.
simplify(true_sol - (4/1.3)*(exp(0.8*t)-exp(-0.5*t))+2*exp(-0.5*t))
ans = I'm not sure what you screwed up in the solve, but your "exact_sol" is not. Not exact, that is. And since I don't see your code, I can only guess what you did wrong on paper. I won't try to do that.
Now let me implement Euler's method. I'll still use your code as a template, but mine will be correct. :)
h = 0.05; % time step
t = 0:h:4; %time range
y = zeros(size(t)); % set y array all 0's to same size as t
y(1) = 2; % set initial condition at time 0 to 2
n = numel(y);
dydt = @(y_i,t_i) 4*exp(0.8*t_i) - 0.5*y_i;
for i=1 : n-1 % for loop to interate through y values for
y(i+1)= y(i)+ h * dydt(y(i),t(i)); % Euler update
end
plot(t,y,'-b') %plot Euler
hold on
fplot(true_sol,[0,4],'r'); % plots the exact solution to this differential equation You can see the two curves as I plotted them are now almost on top of each other. A little larger value of h would have been good perhaps to show off the difference more. But Euler actually did pertty well here.

scott lamb on 26 Nov 2020
Thank you so much for your quick response, i can see my errors now. Though the exact solution was given to me in the question and it says the exact solution and could be determined analytically as : -
(4/1.3)*(exp(0.8*t)-exp(-0.5*t))+2*exp(-0.5*t))
And to write code to plot and compare both the Euler method for y' = 4*exp(0.8*t) - 0.5*y and the exact solution determined analytically?
i have changed
%true_sol = dsolve(diff(y,t,1) == 4*exp(0.8*t) - 0.5*y(t),y(0) == 2);
to
true_sol=(4/1.3)*(exp(0.8*t)-exp(-0.5*t))+2*exp(-0.5*t); %This is the exact solution to dy/dt. determied analytically
and this also seems to work, yeilding the same results? am i totally of the mark or would this now be right given the question? John D'Errico on 26 Nov 2020
syms t
yoursol = (4/1.3)*(exp(0.8*t)-exp(-0.5*t))+2*exp(-0.5*t);
syms y(t)
mysol = dsolve(diff(y,t,1) == 4*exp(0.8*t) - 0.5*y,y(0) == 2);
expand(mysol)
ans = simplify(mysol - yoursol)
ans =
0
So it looks like the two solutions are the same after all after simplification. I probably copied the wrong solution from you.