Trying to use ode45 to solve a second order differential equation with time dependent parameters

Hello,

I am trying to modify some ode45 code to solve second-order differential equation to include time-dependent parameters. So far the code has piece-wise parameters but am interested in how to change the code so that the length parameter is a function of the value of the variable in the second-order differential equation.

Here we set up the second order differential equation

% using ODE45 to solve for theta
ops=odeset('RelTol',1e-10);
% a=1; b=.1; theta_0=pi/2; thetap_0=0; g=9.807; tinit=0; tfinal=98.8
a=2; b=0.6; theta_0=pi/6.04; thetap_0=0; g=9.807; tinit=0; tfinal=3.5
[t,u]=ode45(@exactanswer,[tinit,tfinal],[theta_0,thetap_0],ops,a,b);

Here we set up for ode45:

function uprime=exactanswer(t,u,a,b)
uprime=zeros(2,1);
uprime(1)=u(2);
uprime(2)=-2*lp(t)*u(2)./l(t)-9.807*sin(u(1))./l(t);
end

Here is a part of the piecewise length function

function fval = l(t)
%T = 2*pi*sqrt(a/g);
if (t <= 1.40)
    fval = a;
    T = 2*pi*sqrt(fval/g);
elseif (t > 1.40) && (t <= 1.45)

Is there a way to specify length as a function of time depending on the angle theta?

Thank you.

[Sam Chak]

Hi @Mary Lanzerotti

I'm trying to help and understand the problem. The ODE looks like a pendulum system with a time-varying cord length:

.

You have mentioned that the cord length depends on the time and the angle. Can you mathematically show the functions , and in LaTeX form? It will be easier for me to visualize the geometry of the cord length.

How about describing as well? Because I cannot find a single mathematical function that describe it on this page.

Answer 1

[Walter Roberson]

tinit=0; tfinal=3.5
[t,u]=ode45(@exactanswer,[tinit,tfinal],[theta_0,thetap_0],ops,a,b);

Your time range runs 0 to 3.5 in a single call.

if (t <= 1.40)

You test based on time.

    fval = a;
    T = 2*pi*sqrt(fval/g);

The value you construct is not proportionate to time: the fval is constant up to 1.40 and then suddenly becomes something else.

That means that the derivative of your fval function is discontinuous.

When you use a Runge-Kutta method such as ode45(), the mathematics of the integration require that the calculations are all -- continuous second derivatives. Your code violates that mathematically assumption.

To fix this problem, you have two options:

• change the l(t) function so that the function has continuous second derivatives. For example, interp1() with 'cubic' uses piecewise cubic polynomials that have continuos second derivatives; OR
• each call to ode45() only uses one of the discontinuous branches. You would ask to end integration after 1.40, then invoke ode45() again for 1.40 to 1.45, then use ode45() again as many times as necessary so that there are no discontinuous derivatives within any one ode45() call.

You are asking how to rework as a function of angle. If the angles cannot be calculated directly from time, so that you no longer have one-to-one mappings between times and discontinuities, then you need to introduce an event function into the options. The event function needs to detect that the angle has been reached that signals the boundary of the discontinuity, and the event function has to signal for termination in that case. (Watch out for the possibility that the change could result in the angle regressing temporarily or permanently.)

[Mary Lanzerotti]

Hello, Thank you. Here ode45 is solving for the angle, theta(t). When theta(t) reaches some target value and theta-dot(t) reaches another target value, then the length l(t) will be changing according to some time-dependent function. Since angle is unknown (it is what is being solved for by ode45), how can ode45 be told to solve for angle theta (and theta-dot) using a given length and then use a different length function when theta and theta-dot reach the target values?

[Torsten]

how can ode45 be told to solve for angle theta (and theta-dot) using a given length and then use a different length function when theta and theta-dot reach the target values?

Define an event and interrupt integration when these instances appear.

After the event, restart anew with a modified length function.

Here is a simple example:

https://de.mathworks.com/matlabcentral/answers/1670834-adding-delay-to-a-system-of-differential-equations-using-the-heaviside-function?s_tid=srchtitle

[Mary Lanzerotti]

Hello,

Thank you so much. I am trying to set each call to ode45() to use one of the discontinuous branches. Here I ask to end integration after 0.97, then invoke ode45() again for 0.97 to 1.03. How do we set the outcomes for angle theta and theta-dot as the new inputs [??, ???] in the next call?

a=1; b=.6; theta_0=pi/6; thetap_0=0; g=9.807; tinit=0; tfinal=0.97

[t,u]=ode45(@exactanswer,[tinit,tfinal],[theta_0,thetap_0],ops,a,b);

tinitta = 0.97; tfinaltb = 1.03; tab = (tb+ta)/2;

[t,u]=ode45(@exactanswer,[tinitta,tfinaltb],[??, ???],ops,a,b);

[Mary Lanzerotti]

Torsten,

Thank you. I have been looking at the example and do not understand the example yet.

As a first step, I will try to do the call to ode45() with one of the discontinuous branches first. But I am stuck on how to use the outputs from the first call as inputs into the next call, which is the message above. Thank you kindly.

[Torsten]

If you want to continue integration with the final values of the previous step, use

[t,u]=ode45(@exactanswer,[tinitta,tfinaltb],[u(end,1), u(end,2)],ops,a,b);

But what has changed that made the interruption necessary ?

[Mary Lanzerotti]

Hello Torsten,

Thank you so much. Here we try to do the call to ode45() with time segments. We try to use the outputs from the first call as inputs into the next call. What has changed that made the next time segment necessary is that one of the functions in the differential equation is changing in each time interval. Here we try to do the integration with three time intervals:

function uprime=exactanswer1(t,u,a,b)

uprime=zeros(2,1);

uprime(1)=u(2);

uprime(2)=-2*LP1.*u(2)./L1-9.807*sin(u(1))./L1;

end

function uprime=exactanswer2(t,u,a,b)

uprime=zeros(2,1);

uprime(1)=u(2);

uprime(2)=-2*((-rr).*la.*exp((-(t-ta))*rr)./((1/(la-lb))*la)).*u(2)./((la + la.*exp((-(t-ta))*rr))./((1/(la-lb))*la)+ 0.3333*(la+lb)) - 9.807*sin(u(1))./((la + la.*exp((-(t-ta))*rr))./((1/(la-lb))*la)+ 0.3333*(la+lb));

end

function uprime=exactanswer3(t,u,a,b)

uprime=zeros(2,1);

uprime(1)=u(2);

uprime(2)=-2*LP3.*u(2)./L3-9.807*sin(u(1))./L3;

end

[t1,u1]=ode45(@exactanswer1,[0,ta],[theta_0,thetap_0],ops,a,b);

[t2,u2]=ode45(@exactanswer2,[ta,tb],[u1(end,1),u1(end,2)],ops,a,b);

[t3,u3]=ode45(@exactanswer3,[tb,tc],[u2(end,1),u2(end,2)],ops,a,b);

h=figure

hold on

plot(t1,u1(:,1))

hold on

plot(t2,u2(:,1))

hold on

plot(t3,u3(:,1))

[Torsten]

If you enlarge the list of parameters handed to "exactanswer1/2/3" to compute the uprime values, this should work, shouldn't it ?

[Mary Lanzerotti]

Torsten,

I am sorry, I don't understand. I was wondering if the code segment works, that is, is it calculating correctly the integration in each time segment.

What does 'enlarge the list of parameters handed to 'exactanswer1/2/3' to computer the uprime values' mean?

[Torsten]

I mean how do you plan to get the values of the variables in the functions exactanswer1,2 and 3 that are not defined there, e.g. LP3 and L3 in "exactanswer3" ? Or do you use the nested-function concept so that they are visible from the calling program ?

[Mary Lanzerotti]

Torsten,

Thank you. Here we try to set the variables:

%This is the ODE45 code.

function project_anglefunction

close all

clc

% using ODE45 to solve for theta

ops=odeset('RelTol',1e-10);

a=30; b=.6; theta_0=pi/6; thetap_0=0; g=9.807;

la = 30.0; lb = 0.8*la; lc = 0.6*la; ld = 0.8*la; le = 0.88*la;

lf = 0.70*la; lg = 0.45*la; lh = 0.62*la; li = 0.7*la; lj = 0.4*la; lk = 0.18*la;

ll = 0.3*la; lm = 0.5*la; ln = 0.25*la; lo = 0.2*la

ta = 1.0*sqrt(la); tb = 1.02*sqrt(la); tab = (tb+ta)*sqrt(la)/2;

tc = 1.04*sqrt(la) ; td = 1.06*sqrt(la); tcd = (tc+td)*sqrt(la)/2

te = 1.4*sqrt(la) ; tf = 1.45*sqrt(la); tef = (te+tf)*sqrt(la)/2

rr = 90

function uprime=exactanswer1(t,u,a,b)

uprime=zeros(2,1);

uprime(1)=u(2);

uprime(2)=-2*LP1.*u(2)./L1-9.807*sin(u(1))./L1;

end

function uprime=exactanswer2(t,u,a,b)

uprime=zeros(2,1);

uprime(1)=u(2);

uprime(2)=-2*((-rr).*la.*exp((-(t-ta))*rr)./((1/(la-lb))*la)).*u(2)./((la + la.*exp((-(t-ta))*rr))./((1/(la-lb))*la)+ 0.3333*(la+lb)) - 9.807*sin(u(1))./((la + la.*exp((-(t-ta))*rr))./((1/(la-lb))*la)+ 0.3333*(la+lb));

end

function uprime=exactanswer3(t,u,a,b)

uprime=zeros(2,1);

uprime(1)=u(2);

uprime(2)=-2*LP3.*u(2)./L3-9.807*sin(u(1))./L3;

end

L1 = la

L2 = (la + la.*exp((-(t-ta))*rr))./((1/(la-lb))*la)+ 0.3333*(la+lb)

L3 = lb

L4 = (lb + lb.*exp((-(t-tc))*rr))./((1/(lb-lc))*lb)+ 0.286*(lb+lc)

LP1 = 0

LP2 = (-rr).*la.*exp((-(t-ta))*rr)./((1/(la-lb))*la)

LP3 = 0

LP4 = (-rr).*lb.*exp((-(t-tc))*rr)./((1/(lb-lc))*la)

[t1,u1]=ode45(@exactanswer1,[0,ta],[theta_0,thetap_0],ops,a,b);

[t2,u2]=ode45(@exactanswer2,[ta,tb],[u1(end,1),u1(end,2)],ops,a,b);

[t3,u3]=ode45(@exactanswer3,[tb,tc],[u2(end,1),u2(end,2)],ops,a,b);

h=figure

hold on

plot(t1,u1(:,1))

hold on

plot(t2,u2(:,1))

hold on

plot(t3,u3(:,1))

end

Answer 2

[Torsten]

Results as expected ?

%This is the ODE45 code.

function project_anglefunction
close all
clc  
% using ODE45 to solve for theta
ops=odeset('RelTol',1e-10);
a=30; b=.6; theta_0=pi/6; thetap_0=0; g=9.807; 
la = 30.0; lb = 0.8*la; lc = 0.6*la; ld = 0.8*la; le = 0.88*la; 
lf = 0.70*la; lg = 0.45*la; lh = 0.62*la; li = 0.7*la; lj = 0.4*la; lk = 0.18*la; 
ll = 0.3*la; lm = 0.5*la; ln = 0.25*la; lo = 0.2*la
ta = 1.0*sqrt(la); tb = 1.02*sqrt(la); tab = (tb+ta)*sqrt(la)/2; 
tc = 1.04*sqrt(la) ; td = 1.06*sqrt(la); tcd = (tc+td)*sqrt(la)/2
te = 1.4*sqrt(la) ; tf = 1.45*sqrt(la); tef = (te+tf)*sqrt(la)/2
rr = 90;
L1 = la;
L3 = lb
LP1 = 0
LP3 = 0
[t1,u1]=ode45(@exactanswer1,[0,ta],[theta_0,thetap_0],ops);
[t2,u2]=ode45(@exactanswer2,[ta,tb],[u1(end,1),u1(end,2)],ops);
[t3,u3]=ode45(@exactanswer3,[tb,tc],[u2(end,1),u2(end,2)],ops);
h=figure
hold on
plot(t1,u1(:,1))
hold on
plot(t2,u2(:,1))
hold on
plot(t3,u3(:,1))
    function uprime=exactanswer1(t,u)
        uprime=zeros(2,1);
        uprime(1)=u(2);
        uprime(2)=-2*LP1.*u(2)./L1-9.807*sin(u(1))./L1;
    end

    function uprime=exactanswer2(t,u)
        uprime=zeros(2,1);
        uprime(1)=u(2);
        uprime(2)=-2*((-rr).*la.*exp((-(t-ta))*rr)./((1/(la-lb))*la)).*u(2)./((la + la.*exp((-(t-ta))*rr))./((1/(la-lb))*la)+ 0.3333*(la+lb)) - 9.807*sin(u(1))./((la + la.*exp((-(t-ta))*rr))./((1/(la-lb))*la)+ 0.3333*(la+lb));
    end

    function uprime=exactanswer3(t,u)
        uprime=zeros(2,1);
        uprime(1)=u(2);
        uprime(2)=-2*LP3.*u(2)./L3-9.807*sin(u(1))./L3;
    end

end 

[Mary Lanzerotti]

Torsten,

Good afternoon and thank you kindly. Yes, these results are what the code was producing. Thank you so much!

[Torsten]

My pleasure.

[Sam Chak]

Hi @Mary Lanzerotti

@Torsten's workaround to deal with the piecewise function in solving the ODE is brilliant. When I inspect your piecewise function of the Length, , I discover that there are two discontinuities at and , as circled in the 1st figure below. Although the Reviewer or Examiner may not be able to detect these (unless you submit your code), the discontinuities do not represent the actual physical science.

If there exists a smooth time-dependent function that decribes the variation of the length, as shown in the 4th figure, would you consider using the time-continuous smooth function that is also integrable with ode45? Rest assured that steepness can be adjusted.

Alternatively, consider adding the time-varying in the ODE as well as uprime(3), if you know the governing equation of according to the real physical science.

la = 30.0; 
lb = 0.8*la; 
rr = 90;
ta = 1.0*sqrt(la); 
tb = 1.02*sqrt(la); 
tc = 1.04*sqrt(la); 

t1 = linspace(0, ta);
t2 = linspace(ta, tb);
t3 = linspace(tb, tc);

L1 = la*ones(1, 100);
L2 = ((la + la*exp(rr*(ta - t2)))/(la/(la - lb)) + 0.3333*(la + lb));
L3 = lb*ones(1, 100);

plot(t1, L1, t2, L2, t3, L3)