How to integrate this function using ode23?

The given function is as follows:
I've written the following code:
clc
clear all
tspan = [0 10];
y0 = 0;
%[t,y] = ode23(@(t,y) 2*t, tspan, y0);
Minc = 67;
Mp = 12;
dfe = 7200;
din = 4040;
Dp = 5e-5;
Cl = 0.045;
Cleq = 0.005;
z = (Minc/(100*Mp))*(dfe/din)*Dp*(Cl - Cleq) ;
disp( z)
[t,y] = ode23(@(t,y) (1/y)*(Minc/(100*Mp))*(dfe/din)*Dp*(Cl - Cleq), tspan, y0);
plot(t,y,'-o');

8 comentarios

Ameer Hamza
Ameer Hamza el 14 de Mzo. de 2020
Editada: Ameer Hamza el 14 de Mzo. de 2020
Your code is correct. But the behavior shown by ode23 and other function of ode-suite is quite strange. We know that even an analytical solution exists for your ODE with given initial conditions. But somehow, ode23 is calculating all the values of y to be NaN. After a bit of fiddling, it can be seen that the issue is caused by setting the initial condition to be zero. If you increase the initial condition even by a ridiculously tiny amount, the ode45 will give a solution. For example, change the initial condition to be
y0 = 1e-308; % minimum possible exponent for double
Maybe this is a bug. Maybe someone else here can come up with a valid reason why it should happen.
Anshuman S
Anshuman S el 14 de Mzo. de 2020
It works! btw I'm using ode23, shall I use ode45?
Ameer Hamza
Ameer Hamza el 14 de Mzo. de 2020
Editada: Ameer Hamza el 14 de Mzo. de 2020
For such problems, both ode23 and ode45 should work fine. I mistakenly wrote ode45 in my comment. I wanted to write ode23. I have also added this as an answer so that anyone facing this error in the future can easily find the fix.
What about this equation:
my code:
clc
clear all
tspan = [0 10];
y0 = 1e-308;
Dp = 5e-5;
Cl = 0.045;
Cleq = 0.005;
Cinl = 0.05;
Cin = 0.1;
alpha = 2*((Cinl - Cl)/(Cin - Cinl));
[t,y] = ode23(@(t,y) (-(alpha*Dp/(2*y)) - (alpha/2)*sqrt(Dp/(3.1415*t))), tspan, y0);
plot(t,y,'-o');
I am not sure about this one. Changing the initial condition to some other values helps but I am not sure this is what you want
tspan = [0.5 10];
y0 = 0.5;
Using a stiff solver also seems to help
[t,y] = ode15s(@(t,y) (-(alpha*Dp/(2*y)) - (alpha/2)*sqrt(Dp/(3.1415*t))), tspan, y0);
But In this case, this might be the property of the differential equation itself.
Anshuman S
Anshuman S el 14 de Mzo. de 2020
Oh, my bad!
It was supposed to be a non-zero intial condition! As rightly stated by you, it's the property of the differential equation itself!
Anshuman S
Anshuman S el 14 de Mzo. de 2020
But why tspan = [0.5 10], why not start from 0
tspan = [0 10]
I was trying different values and this is the one I pasted here. This following one works too
tspan = [1e-308 10];

Iniciar sesión para comentar.

 Respuesta aceptada

Ameer Hamza
Ameer Hamza el 14 de Mzo. de 2020
Although I couldn't figure out the actual cause of this issue, the error is caused by setting the initial condition to zero. It can be resolved by giving a tiny value to the initial condition. For example,
y0 = 1e-308; % minimum possible exponent for double

Más respuestas (0)

Categorías

Más información sobre Numerical Integration and Differential Equations en Centro de ayuda y File Exchange.

Productos

Versión

R2019b

Etiquetas

Preguntada:

el 14 de Mzo. de 2020

Comentada:

el 14 de Mzo. de 2020

Community Treasure Hunt

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

Start Hunting!

Translated by