Using Trapezodial rule Matlab

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')

[Jan]

Please explain, what "I am kinda stuck" explicitly means. What is your question?

Answer 1

[James Tursa]

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)