How can I solve this ode problem?

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

[Jan]

" I think that my problem is stiff" - why do you assume this?

" How can I explain this result?" - explain what? Do you expect something else? If so, what?

[Rachele Agnese Cozzetti]

it' hard for ode45 solve this problem, i have a warming that say me that the solver stops at time t=2.1219*10^2 and i have dense solution, this occours about a stiff problem , it's teue?

[Jan]

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]

time: tspan=[0 250]

y0=[1 sqrt(2)]

Answer 1

[Alan Stevens]

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

[Rachele Agnese Cozzetti]

how i can do this?

[Alan Stevens]

It really depends on what physics you want to insert! However, if you assume the acceleration stops when it reduces to zero (which happens when u ~ 54) then you could do the following:

y0 = [1 sqrt(2)];
a=0; b=250; 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( ~,y )

    if y(1) > 54
        y(1) = 54;
    end

    dyode=[y(2);
           (-1/(y(1)^2))+(0.012/((60-y(1))^2))];
end

[Rachele Agnese Cozzetti]

Thank you! Now it's clear!

[Rachele Agnese Cozzetti]

I need your help. I've calculated in this way h=min(diff(t)) to find the discretization step to use with other method cause of I've to do the comparision between method. I calculate th solution but the length is different beetween ode45 that is my reference and other method, ho i can solve this problem?

Thank you

[Alan Stevens]

I'd need to see the coding for your other method.

[Rachele Agnese Cozzetti]

Can I have your email,please?

[Rachele Agnese Cozzetti]

The code is enough long.

[Alan Stevens]

No, just post your problems/coding/data etc. here - you have a much better chance of having your problems solved that way, as more than one person will see them.

[Rachele Agnese Cozzetti]

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

[Rachele Agnese Cozzetti]

%% 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')

%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
%
t0=tspan(1);
tf=tspan(2);
t=t0:h:tf;
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

[Rachele Agnese Cozzetti]

This is one method, than I've used: Euler implicit, Crank-Nicolson, Heun but the length of y is the same. I've to calculate the error between Runge-Kutta and other method, butthe length is not the same.

[Rachele Agnese Cozzetti]

This is my project, see rocket to the moon 1D.

[Alan Stevens]

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]

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

[Alan Stevens]

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