Info

La pregunta está cerrada. Vuélvala a abrir para editarla o responderla.

which differential equation solver should I use?

1 visualización (últimos 30 días)
February
February el 15 de Feb. de 2015
Cerrada: MATLAB Answer Bot el 20 de Ag. de 2021
I have a differential equation of the form following:
dy(t)/dt = ay(t) + bx(t) + cx(t)/dt
y(t) is the output I am looking for, when x(t) is a given input (discrete samples from collected data). a, b, c are constants.
I followed ode45 instructions for differential equation of time dependent variables. To test out the code, I set x(t) = sin(t). Accordingly I set x(t)/dt = cos(t). this probably is not the way I will test out the real data since they are noisy discrete samples. I need to figure out how to find calculate x(t)/dt for those but that's another question.
I got a straight linear line with positive slope for output y(t). I was expected oscillating function like modified sin(t) and cos(t). I tried on wolfram alpha and that displayed oscillating (damping) graph for y(t), so I think I did not solve this equation correctly on MATLAB. I am pretty sure I implemented ode45 correctly, so I am wondering I may be using a wrong differential equation solver.
Could you give me any feedback? Thanks.
--------------------------------------
t = linspace(0,pi*20,pi*20/0.01);
Pa = sin(t)';
dPa = cos(t)';
f = 4 * -2.6021e-07 * ones(size(Pa,1),1);
g = 4 * 2.6021e-07 * Pa + .6 * dPa;
Tspan = [0 pi*20];
IC = 5 * ones(size(Pa,1),1);
t = t';
Tspan = Tspan';
[T Pic] = ode45(@(t,Pic) f1(t,Pic,t,f,t,g),Tspan, IC');
plot(t,Pa,'b') % input
hold on;
plot(T,Pic(:,45)) % one possible output
%
and function f1 is defined as following:
function dydt = f1(t,y,ft,f,gt,g)
dydt = f.*y + g;
end
  3 comentarios
Torsten
Torsten el 17 de Feb. de 2015
You have to supply g at the times indicated by the variable "t" in dydt. Thus your function routine must look like
function dydt = f1(t,y,ft,f,gt,g)
f = 4*(-2.6021e-07);
Pa = sin(t);
dPa = cos(t);
g = 4*2.6021e-07*Pa + .6*dPa;
dydt=f*y+g;
end
And you only have to supply one value as initial condition for y, not a vector of length 5*size(Pa).
Best wishes
Torsten.
February
February el 20 de Feb. de 2015
Thank you, Torsten. I finally figured out.

Respuestas (0)

La pregunta está cerrada.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by