Backward Eurler and Forward Eurler on matrices in matlab
4 次查看(过去 30 天)
显示 更早的评论
I'm trying to to solve the sysmtem of differential equations where I'm given initial condition x(0) = X0 and matrix A and use forward Eurler and Backward Eurler to find X.
So I tried writing the code as follows
load ODE_Data %load the matrix A and initial condition X0
dT = 6/500; %sec
t_be = 0:dT:6;
Xbe = zeros(5,length(t_be));
Xbe(:,1) = X0;
for t=2:length(t_be)
Xbe(:,t) = ((1/dT)-A)\(Xbe(:,t-1)/dT);
end
for Backward Eurler using the equation
and
load ODE_Data %load the matrix A and initial condition X0
%Name your time points as t_fe and Name your solution vector as Xfe.
dT = 6/500; %sec
t_fe = 0:dT:6;
Xfe = zeros(5,length(t_fe));
Xfe(:,1) = X0;
for t=2:length(t_fe)
b =Xfe(:,t-1);
c=dT*(A+(1/dT));
Xfe(:,t) = dT*(A+(1/dT))*Xfe(:,t-1);
end
for Forward Eurler using the equation
but the anwsers for both quickly goes to infinity and I'm really lost on why that is.
0 个评论
回答(1 个)
Michael Hubatka
2022-1-6
First, make sure that
eig(A) < 0
otherwise, the solution should in fact go to infinity.
The real problem in your code is the calculation of
1/dT - A
A + 1/dT
This will modify all entries in A, the correct code is
eye(size(A))./dT - A
A + eye(size(A))./dT
which will modify only the diagonal.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!