How can I solve this ode problem?

Question

Rachele Agnese Cozzetti 2021-3-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/764911-how-can-i-solve-this-ode-problem

评论： 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 个评论
显示 2更早的评论隐藏 2更早的评论

Jan 2021-3-7

编辑：Jan 2021-3-7

There are two poles at u==0 and u==60. There the 2nd derivative gets infinite.

What are the time limits? How did you define y0?

Rachele Agnese Cozzetti 2021-3-7

time: tspan=[0 250]

y0=[1 sqrt(2)]

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2021-3-7

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/764911-how-can-i-solve-this-ode-problem#answer_641186

You get a singularity (divide by zero) when u = 60, giving infinite acceleration. You need to include a guard against this.

15 个评论
显示 13更早的评论隐藏 13更早的评论

Rachele Agnese Cozzetti 2021-3-9

在 MATLAB Online 中打开

%% MODELLO SEMPLIFICATO AD UNA DIMENSIONE:ROCKET TO THE MOON
clc,clear all, close all
%SISTEMA ODE f1 NON LINEARE
a=0; b=250; tspan=[a,b];   %intervallo di integrazione
y0=[1,sqrt(2)];             %condizioni iniziali
%% RUNGE-KUTTA 4,5 
opt = odeset('Stats','on');
tic
[t, y] = ode45(@f1,tspan,y0,opt);
time_RK=toc;
N_RK=length(t); %numeri di passi per calcolare la soluzine
figure,
plot(t,y(:,1)),xlabel('t'),ylabel('y1(t)'), title('Metodo Runge-Kutta 4,5')
h=min(diff(t));
%% RUNGE-KUTTA 2,3
opt = odeset('Stats','on');
tic
[t23, y23] = ode45(@f1,tspan,y0,opt);
time_RK23=toc;
N_RK23=length(t); %numeri di passi per calcolare la soluzine
figure,
plot(t23,y23(:,1)),xlabel('t'),ylabel('y1(t)'), title('Metodo Runge-Kutta 2,3')
%% STIFFNESS
%Ho un problema non lineare quindi per verificare che si tratta di un problema
%stiff vado a calcoloare la matrice Jacobiana e poi calcolo il rapporto di
%stiffness dato dal rapporto tra l'autovalore massimo e l'autovalore minimo
j1=jac1(t,y);
lambda=eig(j1);
minimo=min(abs(lambda));
massimo=max(abs(lambda));
rs=abs(massimo)/abs(minimo); %quoziente di stiffness
%Ho trovato che rs è pari a uno, problema non stiff
%% EULERO ESPLICITO
tic
[t_ee, y_ee] = eul_esp(@f1,tspan,y0,h); 
time_ee=toc;
N_EE=length(t_ee);
figure,
plot(t,y(:,1),t_ee,y_ee(:,1)), xlabel('t'), ylabel('y(t)'), title('Metodo di Eulero esplicito'), legend('RK 45','eulero esplicito')
%% METODO DI HEUN-metodo esplicito
tic
[t_h, y_h] =heun(@f1,tspan,y0,h);
time_heun=toc;
N_H=length(t_h);
figure,
plot(t,y(:,1),t_h,y_h(:,1)), xlabel('t'), ylabel('y(t)'), title('Metodo di Heun'), legend('RK 45','HEUN')
%% Metodo di Crank-Nicolson-implicito
tic
[t_cn, y_cn] = crank_nic(@f1,@jac1,tspan,y0',h);
time_cn=toc;
%numero passi Crank-Nicolson al variare di h
N_CN=length(t_cn) ;
figure,
plot(t,y(:,1),t_cn,y_cn(:,1)), xlabel('t'), ylabel('y(t)'), title('Metodo di Crank-Nicolson'),legend('RK 45','C_N')
%% Eulero implicito
tic
[t_ei, y_ei] = eul_imp(@f1,@jac1,tspan,y0',h);
time_ei=toc;
N_EI=length(t_ei);
figure,
plot(t,y(:,1),t_ei,y_ei(:,1)),xlabel('time'), ylabel('y(t)'), title('Metodo di Eulero implicito'),legend('RK 45','e_imp')
%% In un unico plot faccio il confronto tra le le varie soluzioni
figure,
p=plot(t,y(:,1),'-*m',t23,y23(:,1),'--',t_ee,y_ee(:,1),'--',t_h,y_h(:,1),'--',t_cn,y_cn(:,1),'--',t_ei,y_ei(:,1),'--'), xlabel('time'), ylabel('y(t)'), title('ROCKET TO THE MOON 1D'),legend('RK 45','RK 23','EE','Heun','CN','EI')
%% CALCOLO ERRORE
% e_e=norm(y-y_ee)/norm(y)
% e_h=norm(y-y_h)/norm(y)
% e_ei=norm(y-y_cn)/norm(y)
% e_cn=norm(y-y_ei)/norm(y)
%% Confronto in termini di tempo computazionale
T=table(time_RK,time_RK23, time_ee, time_heun, time_cn, time_ei)

Alan Stevens 2021-3-9

在 MATLAB Online 中打开

Try the following

%% EULERO ESPLICITO
y0 = [1 sqrt(2)];
a=0; b=250;
N = 10000; h = (b - a)/N; % N needs to be larger for an accurate Euler
tspan= a:h:b;
[t, y] = ode45(@f1,tspan,y0,odeset('RelTol',1e-10));
N_RK=length(t);
[t_ee, y_ee] = eul_esp(@f1,tspan,y0,h); 
N_EE=length(t_ee);
figure,
plot(t,y(:,1),'b--',t_ee,y_ee(:,1),'r'),
axis([a b 0 60])
xlabel('t'), ylabel('y(t)'), title('Metodo di Eulero esplicito'),
legend('RK 45','eulero esplicito')
%eul_esp
function [t, y]=eul_esp(funz,tspan,y0,h)
%
% Funzione che calcola la soluzione di un'equazione differenziale 
% con il metodo di Eulero esplicito.
% INPUT:
%		funz	stringa contenente il nome del file in cui
%				inserire l'equazione differenziale.
%		tspan   [t0,tf]	istante di tempo iniziale e finale
%		y0		condizione iniziale
%		h		passo temporale
% OUTPUT:
%		t		vettore contenente un insieme di punti interni
%				all'intervallo desiderato
%		y		vettore contenente la soluzione nei punti del vettore t
%
t=tspan;
n=length(t);
y(1,:)=y0;
for i=1:n-1
   ff= funz(t(i),y(i,:));
   y(i+1,:)=y(i,:)+h*ff';
end
end
function dyode = f1( ~,y )
    if y(1) > 54
        y(1) = 54;
    end
    dyode=[y(2);
           (-1/(y(1)^2))+(0.012/((60-y(1))^2))];
end

NB: you can attach m-files (using the paper-clip symbol) if you have a large file.

Rachele Agnese Cozzetti 2021-3-9

The prof say me that I've to use as h this : h=min(diff(t)) for other method.

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).

请先登录，再进行评论。

How can I solve this ode problem?

4 个评论
显示 2更早的评论隐藏 2更早的评论

采纳的回答

15 个评论
显示 13更早的评论隐藏 13更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

How can I solve this ode problem?

4 个评论 显示 2更早的评论隐藏 2更早的评论

采纳的回答

15 个评论 显示 13更早的评论隐藏 13更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论

15 个评论
显示 13更早的评论隐藏 13更早的评论