Plotting a solution of a differential equation with Sagemath

Question

I need to solve a differential equation. The solution will depend on $t$ and $q$, and I need to define that $q$ piecewise depending on $t$.

var('k,Tmax,Tmin,w,T0,q'); T=function('T')(t); Te=function('Te')(t);

assume(k>0); assume(Tmax>Tmin); Te(t)=(Tmax+Tmin)/2+(Tmax-Tmin)/2*sin(w*t);

Now this is the differential equation solution:

sol=desolve(diff(T(t),t)-q+k*(T(t)-Te(t)),[T,t],[0,T0]);

The solution with $q=0$ for example would be

sol.subs(Tmax=21.6,Tmin=15.2,k=0.024,q=0,T0=15.6,w=pi/12);

but I need that q to be a model for a heater that's on from 8 AM to 22 PM, and off from 22 PM to 8 AM. So I need to define a $q$ function that if $t mod 24$ is between $8$ and $22$ its value is $0.0504$, and $0$ otherwise. Something like this

$$q(t)=\begin{cases}0.0504 \quad\ \ \ \ \ \ \ if \quad t\ mod\ 24 \in[8,22] \\ 0 \ \ \qquad \qquad otherwise\end{cases}$$

I've been trying with the piecewise function but it's not plotting anything, I always get error messages.

Thanks for your time.

Answer 1

You can express $q(t)$ as a sum of differences of step functions. Also, it's more efficient to solve the differential equation numerically. I assume you want to plot the solution for some number of days (which can be specified in the code).

days=3
Tmax=21.6; Tmin=15.2; k=0.024; T0=15.6; w=pi/12
var('t');
Te=(Tmax+Tmin)/2+(Tmax-Tmin)/2*sin(w*t)

We define $q(t)$:

q = 0.0504*sum(unit_step(t-8-d*24) - unit_step(t-22-d*24) for d in range(days))
plot(q,t,0,days*24)

We define the ODE (also the one with $q(t)=0$, to compare):

T=function('T')(t);
ode0 = diff(T,t) == -k*(T-Te)
ode1 = diff(T,t) == q-k*(T-Te)

Finally we solve and plot:

sol0=desolve_rk4(ode0, T, ivar=t, ics=[0,T0], step=0.1, end_points=[0,days*24], output='plot', xmin=0,xmax=days*24)
sol1=desolve_rk4(ode1, T, ivar=t, ics=[0,T0], step=0.1, end_points=[0,days*24], output='plot', xmin=0,xmax=days*24, color='red')
sol0 + sol1