using ode45

Question

1 个投票

Dear all,

I would like to solve the following system of differential equations with input with Matlab. I tried with ODE45. Unfortunately, I always get a very strange behaviour of the solution, almost insensitive to the input. The values of the parameters (taus tauf etc..) are taken from a work where integration of the same equations was succesfull, and are supposed to be meaningful. The input should be an impulse at given time. The solution should look null before the impulse, increasing afterwards, reaching zero again after about 20 sec.

It is probably some stupid mistake but I cannot get it right. Any suggestion?

Thanks in advance, Sara

%%---------------------------------

pnts = 7680;

dt = 0.0039;

x0 = zeros(1,2);

u = zeros(pnts,1);

u(768) = 1;

t = [1 pnts]*dt;

[T, SOL] = ode45(@fx, t, x0,[],u,dt);

%% ----------------------------

function dx = fx(t,x,u,dt)

ui = u(ceil(t/dt));

x = x';

taus = 0.8;

tauf = 0.4;

epsilon = 0.5;

dx(:,1) = epsilon*ui-x(:,1)./taus-(x(:,2)-1)./tauf;

dx(:,2) = x(:,2);

dx = dx';

end

%%---------------------------------------------------------------

6 个评论
显示 4更早的评论隐藏 4更早的评论

Andrew Newell 2012-1-24

What is u? You have it equal to zero except at one point.

Sara 2012-1-24

u is an impulse that is one only at one time point and zero otherwise. I did try to use eps instead of zero, but with no difference in the results. I also try with u=randn(pnts,1), but still i do not get the output I expect. I am wondering whther it is a problem of transient behaviour of the equations or ode45 is not the appropriate solver to use...

请先登录，再进行评论。

请先登录，再回答此问题。

Follow Question

Answer 1

Andrew Newell 2012-1-24

在 MATLAB Online 中打开

1 个投票

See How Do I Obtain an Impulse Response? for the right way to model impulses.

EDIT: Here is how to apply it to your problem, assuming that your impulse occurs at t=3.

First, integrate the impulse-free equation for the first three seconds. Your equations can be represented by a simple online function:

taus = 0.8;
tauf = 0.4;
fx = @(t,x) [-x(1)./taus-(x(2)-1)./tauf; x(1)];
x0 = zeros(1,2);
tspan = [0 3];
[T, SOL] = ode45(fx, tspan, x0);

Now add the impulse at t=3 and use the result as the initial condition for the next seven seconds:

epsilon = 0.5;
x0 = SOL(end,:) + [epsilon 0];
tspan = [3 10];
[T1, SOL1] = ode45(fx,tspan,x0);

Finally, combine the two segments to get your solution:

T = [T; T1]; SOL = [SOL; SOL1];

4 个评论
显示 2更早的评论隐藏 2更早的评论

Sara 2012-1-25

I ran the code you suggested a couple of times. IN this figure (http://fileslap.com/93x/ComparisonImpulseTime) you can see the results (green is x1, blue is x2; thick is with an impulse at 3; dashed is with impulse at 10). The dashed lines look correct. So, my understanding is that the system as some kind of transient, that I have to wait before giving the impulse. I guess I can calculate this transient time from the time costants, but this is not matlab-related I would say.

Nevertheless, I do have an additional question: if I have to simulate a train of impulses as input do I have to break the calculation into segments for each input? Is by any chance the "event" feature of the ode45 related to this matter?

Thanks

Sara

Andrew Newell 2012-1-25

The only difference between the dashed output and that of my code is the time intervals. Just replace 3 by 10 for the impulse time and 10 by 20 for the total time. As for a train of impulses, if you know the time of each impulse you can create a loop to piece together multiple time segments.

请先登录，再进行评论。

using ode45

6 个评论
显示 4更早的评论隐藏 4更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

类别

标签

Community Treasure Hunt

using ode45

6 个评论 显示 4更早的评论 隐藏 4更早的评论

采纳的回答

4 个评论 显示 2更早的评论 隐藏 2更早的评论

更多回答（0 个）

类别

标签

另请参阅

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论