How can I solve this ode problem?
2 ビュー (過去 30 日間)
古いコメントを表示
Rachele Agnese Cozzetti
2021 年 3 月 6 日
コメント済み: Alan Stevens
2021 年 3 月 9 日
I've a probem of rocket to the moon in one dimension, it's described by a non linear 2nd order ODE with two initial condition:
a=0; b=1000; tspan=[a b];
[t, y] = ode45(@f1,tspan,y0,odeset('RelTol',1e-10));
N_RK=length(t)
plot(t,y(:,1),'or')
xlabel('t')
ylabel('x(t)')
title('ROCKET TO THE MOON 1D, ode45')
%f1 is my ode
function dyode = f1( t,y )
dyode=[y(2);
(-1/(y(1)^2))+(0.012/((60-y(1))^2))];
end
This is my plot and I think that my problem is stiff but the ratio of stifness is 1. How can I explain this result?
I've found ratio stiffness in this way:
j1=jac1(y);
lambda=eig(j1);
minimo=min(abs(lambda));
massimo=max(abs(lambda));
rs=abs(massimo)/abs(minimo) %quoziente di stiffness
if rs>1
disp('Il problema è stiffness')
else
disp('Il problema non è stiffness')
end
%jac1 Jacobian matrix
function dyy= jac1( y )
dyy=[0 1;
2/(y(1)^3)-0.024/((y(1)-60)^3) 0];
end
%Any error in my code? Thank you
4 件のコメント
採用された回答
Alan Stevens
2021 年 3 月 7 日
You get a singularity (divide by zero) when u = 60, giving infinite acceleration. You need to include a guard against this.
15 件のコメント
Alan Stevens
2021 年 3 月 9 日
Ok, but this pretty well guarantees that the length of t and t_ee will be different, as ode45 chooses it's own timestep length. The way to avoid this is to choose a tspan for ode45 that specifies that it outputs times at a constant interval (as in my code above) in which case min(diff(t)) will give the same value.
I don't see why you need to do this though, as, even if the number of timesteps are different between ode45 and Euler, you can still compare the values at a given time by using interpolation (see interp1).
その他の回答 (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!