Using Trapezodial rule Matlab
3 次查看(过去 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!