Solving a system of ODEs whose coefficients are piecewise functions
2 views (last 30 days)
Show older comments
I try to plot the solution of a system of ODE, on [-10,10], for the initial data [0.001 0.001], using the function:
function dwdt=systode(t,w)
if 0< t<1
f = t*(3-2*t);
if -1<t< 0
f=t*(3+2*t);
else
f = 1/t;
end;
if 0< t <1
h=4*t^4-12*t^3+9*t^2-4*t+3;
if -1< t < 0
h=4*t^4+12*t^3+9*t^2+4*t+3;
else
h=0;
end;
beta=0.5+exp(-abs(t));
dwdt=zeros(2,1);
dwdt(1)=-f*w(1)+w(2);
dwdt(2)=-beta*w(1)-f*w(2)+h*w(1)-f*w(1)^2;
end
The coefficients f(t) and g(t) are piecewise functions as follows.
With the commands
tspan = [-10 10];
z0=[0.001 0.001];
[t,z] = ode45(@(t,z) systode(t,z), tspan, z0);
figure
plot(t,z(:,1),'r');
I get the message
tspan = [-10 10];
↑
Error: Invalid use of operator.
Where could be the mistake? I am also not sure that I defined correctly the functions f, h, β.
0 Comments
Accepted Answer
Torsten
on 28 Dec 2021
function main
T = [];
Z = [];
z0=[0.001 0.001];
tspan1 = [-10 -1];
iflag = 1;
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan1, z0);
T = vertcat(T,t);
Z = vertcat(Z,z);
tspan2 = [-1 0];
iflag = 2;
z0 = [z(end,1) z(end,2)];
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan2, z0);
T = vertcat(T,t);
Z = vertcat(Z,z);
tspan3 = [0 1];
iflag = 3;
z0 = [z(end,1) z(end,2)];
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan3, z0);
T = vertcat(T,t);
Z = vertcat(Z,z);
tspan4 = [1 10];
iflag = 4;
z0 = [z(end,1) z(end,2)];
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan4, z0);
T = vertcat(T,t);
Z = vertcat(Z,z);
figure
plot(T,Z(:,1),'r');
end
function dwdt = systode(t,w,iflag)
if iflag == 1
f = 1/t;
h = 0;
beta=0.5+exp(t);
elseif iflag == 2
f = t*(3+2*t);
h = 4*t^4+12*t^3+9*t^2+4*t+3;
beta=0.5+exp(t);
elseif iflag == 3
f = t*(3-2*t);
h = 4*t^4-12*t^3+9*t^2-4*t+3;
beta=0.5+exp(-t);
elseif iflag == 4
f = 1/t;
h = 0;
beta = 0.5+exp(-t);
end
dwdt=zeros(2,1);
dwdt(1)=-f*w(1)+w(2);
dwdt(2)=-beta*w(1)-f*w(2)+h*w(1)-f*w(1)^2;
end
3 Comments
Torsten
on 28 Dec 2021
Put the code in a file, name it main.m and load it into MATLAB. Run it and the first function will be plotted in the interval [-10:10]. The command is
figure
plot(T,Z(:,1),'r');
More Answers (1)
Cris19
on 29 Dec 2021
4 Comments
Torsten
on 29 Dec 2021
Edited: Torsten
on 29 Dec 2021
This is the more usual formulation for your ODE system.
Seems it was not absolutely necessary to interrupt integration at the points where the piecewise functions were glueed together.
function main
tspan = [-10 10];
z0 = [0.001 0.001];
[t,z] = ode15s(@(t,z) systode(t,z), tspan, z0);
for i=1:numel(t)
dz(:,i) = systode(t(i),z(i,:));
end
figure
plot(t,z(:,1),'r')
figure
plot(t,dz(1,:).','g')
end
function dwdt = systode(t,w)
if abs(t) < 1
f = t*(3-2*abs(t));
else
f = 1/t;
end
if t>=0 && t < 1
h=4*t^4-12*t^3+9*t^2-4*t+3;
elseif t>-1 && t <= 0
h=4*t^4+12*t^3+9*t^2+4*t+3;
else
h=0;
end
beta=0.5+exp(-abs(t));
dwdt=zeros(2,1);
dwdt(1)=-f*w(1)+w(2);
dwdt(2)=-beta*w(1)-f*w(2)+h*w(1)-f*w(1)^2;
end
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!