Using Trapezodial rule Matlab

3 次查看(过去 30 天)
Abdul Qadoos
Abdul Qadoos 2021-11-14
编辑: Jan 2022-4-11
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 个评论
Jan
Jan 2021-11-14
编辑:Jan 2022-4-11
Please explain, what "I am kinda stuck" explicitly means. What is your question?

请先登录,再进行评论。

回答(1 个)

James Tursa
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)

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by