ode45 output incorrect - possible incorrect function definition?

1 visualización (últimos 30 días)
Jess
Jess el 25 de Mayo de 2014
Respondida: Star Strider el 25 de Mayo de 2014
I have two functions, one to define the differential equation and one to perform ode45 and solve the equation numerically. However, ode45 seems to give a very suspicious output which after solving the equation by hand I can see clearly that it is incorrect. I feel that I have somehow defined the equation incorrect but cannot see how. If possible could someone just have a once over the code to see if I am missing something anywhere, thanks.
Function to define the differential equation:
function dndt = productdiffeqn(t,n,par)
a = par(1);
m = par(2);
c = par(3);
dndt = ((a*10*(c-n))/(m+c-n));
end
Function to solve it:
function ndata = productsol(par, tyear)
n0 = 0;
[tdata,ndata] = ode45(@productdiffeqn,[tyear],n0,[],[par]);
tdata
end
Called with:
productsol([1,1,1],[tyear])
tdata =
0
1
2
3
4
6
12
ans =
0
1.0007
0.9995
1.0000
0.9991
0.9994
0.9999
  2 comentarios
Star Strider
Star Strider el 25 de Mayo de 2014
What differential equation did you start with?
We need to know that and other details of what you’re doing before we can figure out what’s wrong.
Jess
Jess el 25 de Mayo de 2014
Editada: Jess el 25 de Mayo de 2014
dn/dt = (10a(c-n))/(m+c-n)
Constants a,c,m are all >0. Initial value of n: n(0) = 0.

Iniciar sesión para comentar.

Respuestas (1)

Star Strider
Star Strider el 25 de Mayo de 2014
The constants a, c, and m all have to be in your workspace, but you can solve your ODE with an anonymous function:
c = 1;
m = 5;
a = 7;
productdiffeqn = @(t,n) (10.*a.*(c-n))./(m+c-n);
tyear = [0 1 2 3 4 6 12];
n0 = 0;
[t, n] = ode45(productdiffeqn, tyear, n0);
figure(1)
plot(t, n)
grid
I don’t know what your constants are, but this looks similar to what you posted.

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by