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 commentaires

Ameer Hamza
Ameer Hamza le 14 Mar 2020
Modifié(e) : Ameer Hamza le 14 Mar 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 le 14 Mar 2020
It works! btw I'm using ode23, shall I use ode45?
Ameer Hamza
Ameer Hamza le 14 Mar 2020
Modifié(e) : Ameer Hamza le 14 Mar 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 le 14 Mar 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 le 14 Mar 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];

Connectez-vous pour commenter.

 Réponse acceptée

Ameer Hamza
Ameer Hamza le 14 Mar 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

Plus de réponses (0)

Catégories

En savoir plus sur Numerical Integration and Differential Equations dans Centre d'aide et File Exchange

Produits

Version

R2019b

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by