Solving a system of ODEs whose coefficients are piecewise functions
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Akzeptierte Antwort
Torsten
am 28 Dez. 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 Kommentare
Torsten
am 28 Dez. 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');
Weitere Antworten (1)
Cris19
am 29 Dez. 2021
4 Kommentare
Torsten
am 29 Dez. 2021
Bearbeitet: Torsten
am 29 Dez. 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
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!