Solving a system of ODEs whose coefficients are piecewise functions
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Cris19
am 28 Dez. 2021
Beantwortet: Sam Chak
am 21 Jun. 2024
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/846240/image.png)
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 (2)
Sam Chak
am 21 Jun. 2024
Hi @Dehua Kang
The red curve in your image is
and the green curve is
. Below is another demo, with the data from the piecewise smooth functions
and
can be easily obtained from the spreadsheet (Google Sheets or MS Excel).
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1719676/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1719681/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1719686/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1719691/image.png)
%% Piecewise smooth functions (in data form)
tpw = linspace(-10, 10, 201);
fpw = [-1/10, -10/99, -5/49, -10/97, -5/48, -2/19, -5/47, -10/93, -5/46, -10/91, -1/9, -10/89, -5/44, -10/87, -5/43, -2/17, -5/42, -10/83, -5/41, -10/81, -1/8, -10/79, -5/39, -10/77, -5/38, -2/15, -5/37, -10/73, -5/36, -10/71, -1/7, -10/69, -5/34, -10/67, -5/33, -2/13, -5/32, -10/63, -5/31, -10/61, -1/6, -10/59, -5/29, -10/57, -5/28, -2/11, -5/27, -10/53, -5/26, -10/51, -1/5, -10/49, -5/24, -10/47, -5/23, -2/9, -5/22, -10/43, -5/21, -10/41, -1/4, -10/39, -5/19, -10/37, -5/18, -2/7, -5/17, -10/33, -5/16, -10/31, -1/3, -10/29, -5/14, -10/27, -5/13, -2/5, -5/12, -10/23, -5/11, -10/21, -1/2, -10/19, -5/9, -10/17, -5/8, -2/3, -5/7, -10/13, -5/6, -10/11, -1, -27/25, -28/25, -28/25, -27/25, -1, -22/25, -18/25, -13/25, -7/25, 0, 7/25, 13/25, 18/25, 22/25, 1, 27/25, 28/25, 28/25, 27/25, 1, 10/11, 5/6, 10/13, 5/7, 2/3, 5/8, 10/17, 5/9, 10/19, 1/2, 10/21, 5/11, 10/23, 5/12, 2/5, 5/13, 10/27, 5/14, 10/29, 1/3, 10/31, 5/16, 10/33, 5/17, 2/7, 5/18, 10/37, 5/19, 10/39, 1/4, 10/41, 5/21, 10/43, 5/22, 2/9, 5/23, 10/47, 5/24, 10/49, 1/5, 10/51, 5/26, 10/53, 5/27, 2/11, 5/28, 10/57, 5/29, 10/59, 1/6, 10/61, 5/31, 10/63, 5/32, 2/13, 5/33, 10/67, 5/34, 10/69, 1/7, 10/71, 5/36, 10/73, 5/37, 2/15, 5/38, 10/77, 5/39, 10/79, 1/8, 10/81, 5/41, 10/83, 5/42, 2/17, 5/43, 10/87, 5/44, 10/89, 1/9, 10/91, 5/46, 10/93, 5/47, 2/19, 5/48, 10/97, 5/49, 10/99, 1/10];
hpw = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 354/625, 659/625, 909/625, 1104/625, 2, 1359/625, 1449/625, 1544/625, 1674/625, 3, 1674/625, 1544/625, 1449/625, 1359/625, 2, 1104/625, 909/625, 659/625, 354/625, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
plot(tpw, fpw), grid on, xlabel('t'), title('Piecewise function, f(t)')
plot(tpw, hpw), grid on, xlabel('t'), title('Piecewise function, h(t)')
%% Solving the ODEs
tspan = [-10 10];
w0 = [0.001 0.001];
sol = ode45(@(t, w) ode(t, w, tpw, fpw, hpw), tspan, w0);
t = linspace(-10, 10, 2001);
[w, wp] = deval(sol, t); % wp is w-prime (w') = dw/dt
plot(t, w(1,:)), grid on, xlabel('t'), title('Solution, w_{1}(t)')
plot(t, wp(1,:)), grid on, xlabel('t'), title('Derivative, dw_{1}/dt')
%% ODEs
function dwdt = ode(t, w, tpw, fpw, hpw)
f = interp1(tpw, fpw, t); % interpolated piecewise function f(t)
h = interp1(tpw, hpw, t); % interpolated piecewise function h(t)
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
0 Kommentare
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!