Using Trapezodial rule Matlab
5 views (last 30 days)
Show older comments
I have an ode of dy/dt = -y, y(0) = 1. I have been trying to implement the trapezoidal rule to print out the error and error orders but I am kinda stuck. This is what I got so far.
clc
clear all
close all
format long
hold on;
for h=[0.5,0.5/2,0.5/4,0.5/8,0.5/16]
x=0:h:1;
y=[1];
for i=2:length(x)
%y(i)=(1-h)*y(i-1); % FEM
%y(i)=y(i-1)/(1+h); % BEM % just change code when wanting FEM or BEM
y(i)=trapz(exp(-1)); % Method for Trapezodial
end
plot(x,y);
fprintf('h = %f -> error %f -> order %f\n',h,abs(exp(-1)-y(end)),abs(exp(-1)+ h));
end
legend('h=0.5','h=0.5/2','h=0.5/4','h=0.5/8','h=0.5/16')
1 Comment
Answers (1)
James Tursa
on 14 Nov 2021
Edited: James Tursa
on 16 Nov 2021
The trapezoid rule is for situations where you know the function values in advance, but you don't have this situation. I.e., if the right hand side was a function of t only, then you could calculate all the right hand side function values in advance and use the trapezoid rule. But you don't have this situation. Your right hand side is a function of y, which you don't know. What is the actual wording of your assignment?
*** EDIT ***
Upon looking at your ODE a bit closer I see there is a workaround for your particular situation. Consider the following integration using a trapezoid:
y(i+1) = y(i) + h * (ydot(i) + ydot(i+1))/2
For your particular ODE, we can substitute for the ydot values:
y(i+1) = y(i) + h * (-y(i) - y(i+1))/2
Then solve this for y(i+1) in terms of y(i):
y(i+1) + y(i+1) * h/2 = y(i) - y(i) * h/2
y(i+1) * (1 + h/2) = y(i) * (1 - h/2)
y(i+1) = y(i) * (1 - h/2) / (1 + h/2)
So you can code this up to take your time steps. Note that this technique only works because your particular ODE allowed us to solve for y(i+1) easily.
Or you can change the indexing to this of course:
y(i) = y(i-1) * (1 - h/2) / (1 + h/2)
0 Comments
See Also
Categories
Find more on Ordinary Differential Equations 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!