Euler's Method and exact solution plot
101 views (last 30 days)
Show older comments
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
legend('Euler','Exact'); % adds a legend
grid on
0 Comments
Accepted Answer
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)
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))
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.
2 Comments
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)
simplify(mysol - yoursol)
So it looks like the two solutions are the same after all after simplification. I probably copied the wrong solution from you.
More Answers (0)
See Also
Categories
Find more on Function Creation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!