フィルターのクリア

Solving a system of ODEs whose coefficients are piecewise functions

2 ビュー (過去 30 日間)
Cris19
Cris19 2021 年 12 月 28 日
回答済み: Sam Chak 2024 年 6 月 21 日
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, β.

採用された回答

Torsten
Torsten 2021 年 12 月 28 日
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 件のコメント
Torsten
Torsten 2021 年 12 月 28 日
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');
Cris19
Cris19 2021 年 12 月 28 日
Thank you a lot! It works very well.

サインインしてコメントする。

その他の回答 (2 件)

Cris19
Cris19 2021 年 12 月 29 日
One more question, if possible...
I would also need to plot on the same graph the derivative of first component of the solution, dwdt(1). From a previous question posted on this forum, I learned that in the case when the coefficients are not piecewise defined, the code can be:
tspan = [-10 10];
z0=[0.001 0.001];
[t,z] = ode45(@(t,z) systode(t,z), tspan, z0);
for k = 1:numel(t)
dwdt(:,k) = systode(t(k),z(k,:));
end
figure
plot(t,z(:,1),'b')
hold on
plot(t, dwdt(1,:), 'k')
hold off
legend('$x(t)$','$\dot{x}(t)$', 'Interpreter','latex', 'Location','best');
But in the present case, I can not handle it. Would you be kind to help me?
  5 件のコメント
Cris19
Cris19 2021 年 12 月 29 日
編集済み: Cris19 2021 年 12 月 29 日
Yes, indeed, it seems to be easier this way. I think it is possible, since the functions f, h, β are continuous (and also differentiable) on the whole [-10,10].
Dehua Kang
Dehua Kang 2024 年 6 月 21 日
The result of Cris19 is
and the result of Torsten is as follows
so there is something different between the two results. The second program can make my matlab (2023b) no respond.

サインインしてコメントする。


Sam Chak
Sam Chak 2024 年 6 月 21 日
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).
%% 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

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

タグ

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by