Euler's Method and exact solution plot

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

legend('Euler','Exact'); % adds a legend

grid on

Answer 1

[John D'Errico]

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]

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]

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.