Function with form, that changes in time

Hello,

I want to solve a differential equation with ode45, but one of my parameters is a function of time. It changes its form depending on the parameter t; for example, for t>0 and t<1 it's a linear function, but for t>=1 it's constant. How can I programme such parameter in MatLab?

Answer 1

[Jan]

Remember that ode45 is designed to integrate only functions, which can be differentiated smoothly. If the function has discontinuities at t==1, integrate in steps:

1. from 0 to 1 with the first definition using the original initial values
2. from 1 to end with the second one using the final value of step 1 as new initial value

See also: http://www.mathworks.com/matlabcentral/answers/59582#answer_72047

Answer 2

[Michal Mirowski]

Ok, great. Now I've got more complicated one: I need to consider this Mo parameter in an equation:

dy =[(Ut-Rt*y(1)-Cc*y(2))/Lt;
     (Cm*y(1)-Mo)/Jm];

Is there any simple way of doing this? Attached PNG file.

[Walter Roberson]

You need to break that up into sections each place the function changes, if you want to use ode45.

[Michal Mirowski]

Ok, thank You very much :)

Answer 3

[Michal Mirowski]

I wrote:

 function [] = model()
% TECHNOLOGIE INFORMACYJNE - PROJEKT

% OKREŚLENIE WARTOŚCI WSPÓŁCZYNNIKÓW MASZYNY DC
global Rt Lt Cc Cm Jm Ut Mo
Rt = 0.6 ;
Lt = 0.012;
Cc = 1.8;
Cm = 1;
Jm = 1;
Ut = 240;

% IMPLEMENTACJA FUNKCJI ODE45
% CZĘŚĆ PIERWSZA - ROZWIĄZANIE RR NIEZALEŻNEGO OD Mo=f(t)
hold on
[T,Y]=ode45(@rownania,[0 20],[0 0]);
plot(T,Y(:,1),'g-');

% CZĘŚĆ DRUGA - ROZWIĄZANIE RR ZALEŻNEGO OD Mo=f(t)
for t=0:0.1:5
    Mo = 5*t + 5;
    [T,Y]=ode45(@rownania,[0 5],[0 0]);
    plot(T,Y(:,2),'b-.');
end
for t=5:0.1:10
    Mo = 30;
    [T,Y]=ode45(@rownania,[5 10],[5 127.5]);
    plot(T,Y(:,2),'b-.');
end
for t=10:0.1:15
    Mo = -3*t + 60;
    [T,Y]=ode45(@rownania,[10 15],[5 123.3]);
    plot(T,Y(:,2),'b-.');
end
for t=15:0.1:20
    Mo = 15;
    [T,Y]=ode45(@rownania,[15 20],[5 126]);
    plot(T,Y(:,2),'b-.');
end
grid on;
legend('Natężenie prądu twornika','Prędkość kątowa twornika');
end

% RÓWNANIA OPISUJĄCE MODEL MATEMATYCZNY MASZYNY DC

function dy = rownania(~,y)
global Rt Lt Cc Cm Jm Ut Mo
dy =[(Ut-Rt*y(1)-Cc*y(2))/Lt;
     (Cm*y(1)-Mo)/Jm];
end

But it doesn't work properly; can You explain to me, where is the problem? I can't see it :/

Edit: it contains only the last part of the solution (15 to 20), and I don't know how to find the "initiation values" for the 2nd, 3rd and 4th step.

[Walter Roberson]

function dy = rownania(t,y,Rt,Lt,Cc,Cm,Jm,Ut,Mo)
dy =[(Ut-Rt*y(1)-Cc*y(2))/Lt;
     (Cm*y(1)-Mo(t)/Jm];
end

together with

    Mo = @(t) 5*t + 5;
    [T, Y] = ode45(@(t,y)rownania(t,y,Rt,Lt,Cc,Cm,Jm,Ut,Mo),[0 5],[0 0]);
    plot(T,Y(:,2),'b-.');
    yout = Y(end,:);

    Mo = @(t) 30*ones(size(t));
    [T, Y] = ode45(@(t,y)rownania(t,y,Rt,Lt,Cc,Cm,Jm,Ut,Mo),[5 10], yout);
    plot(T,Y(:,2),'b-.');
    yout = Y(end,:);

    Mo = @(t) -3*t + 60;
    [T, Y] = ode45(@(t,y)rownania(t,y,Rt,Lt,Cc,Cm,Jm,Ut,Mo),[10 20], yout);
    plot(T,Y(:,2),'b-.');
    yout = Y(end,:);

and so on.

Please note that for your first call

% CZĘŚĆ PIERWSZA - ROZWIĄZANIE RR NIEZALEŻNEGO OD Mo=f(t)
hold on
[T,Y]=ode45(@rownania,[0 20],[0 0]);
plot(T,Y(:,1),'g-');

that Mo has not been assigned a value yet and so because it is a global would be [] which would give you [] for the second dy result which would have failed. I have not shown how to replace that case because I do not know what is desired there. Perhaps

Mo = @(t) zeros(size(t));

[Michal Mirowski]

With these modifications:

% CZĘŚĆ DRUGA - ROZWIĄZANIE RR ZALEŻNEGO OD Mo=f(t)
for t=0:0.1:5
    Mo = @(t) 5*t + 5;
    [T, Y]=ode45(@(t,y)rownania(t,y,Rt,Lt,Cc,Cm,Jm,Ut,Mo),[0 5],[0 0]);
    plot(T,Y(:,2),'b-.');
    yout = Y(end,:);
end
for t=5:0.1:10
    Mo = @(t) 30*ones(size(t));
    [T, Y] = ode45(@(t,y)rownania(t,y,Rt,Lt,Cc,Cm,Jm,Ut,Mo),[5 10],yout);
    plot(T,Y(:,2),'b-.');
    yout = Y(end,:);
end
for t=10:0.1:15
    Mo = @(t) -3*t + 60;
    [T, Y] = ode45(@(t,y)rownania(t,y,Rt,Lt,Cc,Cm,Jm,Ut,Mo),[10 15],yout);
    plot(T,Y(:,2),'b-.');
    yout = Y(end,:);
end
for t=15:0.1:20
    Mo = @(t) 15*ones(size(t));
    [T, Y] = ode45(@(t,y)rownania(t,y,Rt,Lt,Cc,Cm,Jm,Ut,Mo),[15 20],yout);
    plot(T,Y(:,2),'b-.');
end
grid on;
legend('Natężenie prądu twornika','Prędkość kątowa twornika');
end

% RÓWNANIA OPISUJĄCE MODEL MATEMATYCZNY MASZYNY DC

function dy = rownania(t,y,Rt,Lt,Cc,Cm,Jm,Ut,Mo)
dy =[(Ut-Rt*y(1)-Cc*y(2))/Lt;
     (Cm*y(1)- Mo)/Jm];
end

MatLab doesn't seem to like line:

[T, Y]=ode45(@(t,y)rownania(t,y,Rt,Lt,Cc,Cm,Jm,Ut,Mo),[0 5],[0 0]);

[Michal Mirowski]

I will just add, that I am really grateful for what You're doing and I don't know how to thank You for all Your help!

[Walter Roberson]

Do not loop on t. Those calls I showed are instead of looping.

I missed a ) when I was typing.

function dy = rownania(t,y,Rt,Lt,Cc,Cm,Jm,Ut,Mo)
dy =[(Ut-Rt*y(1)-Cc*y(2))/Lt;
     (Cm*y(1)-Mo(t))/Jm];
end

[Michal Mirowski]

Ok, now I see it.

I also didn't mark, that Mo is a function of t: Mo(t).

Now the whole programme works properly. Thank You very much (again!) for all Your help! I believe I've learned something :)