Info

This question is closed. Reopen it to edit or answer.

which differential equation solver should I use?

1 view (last 30 days)
February
February on 15 Feb 2015
Closed: MATLAB Answer Bot on 20 Aug 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 Comments
Torsten
Torsten on 17 Feb 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 on 20 Feb 2015
Thank you, Torsten. I finally figured out.

Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!