Using Trapezodial rule Matlab
조회 수: 5 (최근 30일)
이전 댓글 표시
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
답변 (1개)
James Tursa
2021년 11월 14일
편집: James Tursa
2021년 11월 16일
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
참고 항목
카테고리
Help Center 및 File Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!