Solve a system of Differential Equations with a piecewise function
29 次查看(过去 30 天)
显示 更早的评论
I'm trying to solve a system of 2 differential equations (with second, first and zero order derivatives)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/190564/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/190565/image.jpeg)
in which there is a piecewise function
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/190566/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/190567/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/190568/image.jpeg)
This problem comes from the analysis of a vibrating system.
The unknowns of the system and the piecewise function are functions of time.
The unknowns are: 1. beta=beta(t) ; 2. x=x(t)
beta and x with one dot at the top are first order derivatives (respect to time).
beta and x with two dots at the top are second order derivatives (respect to time).
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/190569/image.jpeg)
I have to obtain a plot for each of the two unknows and their second and first order derivatives (respect to time).
To solve the problem, I've tried every system: ODE (returned error), dsolve (impossible to find solution), ecc. … without results.
The only way I think it can solve the problem is to separate the various components of the piecewise function and use dsolve to obtain distinct solutions. Unfortunately, in this way code runs for at least 12 hours to plot data.
There is a better way to solve my problem?
6 个评论
Torsten
2018-5-17
if (t>=0)&&(t<=T1)
alpha=alphamax*(10*(t/T1)^3-15*(t/T1)^4+6*(t/T1)^5);
Dalpha=alphamax*(10*3*(t/T1)^2/T1-15*4*(t/T1)^3/T1+6*5*(t/T1)^4/T1;
D2alpha=alphamax*(10*3*2*(t/T1)/T1/T1-15*4*3*(t/T1)^2/T1/T1+6*5*4*(T/T1)^3/T1/T1;
elseif (t>T1)&&(t<=T2)
alpha=alphamax;
Dalpha=0.0;
D2alpha=0.0;
elseif (t>T2)&&(t<=T3)
alpha=alphamax*(10*((T3-t)/T1)^3-15*((T3-t)/T1)^4+6*((T3-t)/T1)^5);
Dalpha=alphamax*(10*3*((T3-t)/T1)^2*(-1/T1)-15*4*((T3-t)/T1)^3*(-1/T1)+6*5*((T3-t)/T1)^4*(-1/T1));
D2alpha=alphamax*(10*3*2*((T3-t)/T1)*(-1/T1)*(-1/T1)-15*4*3*((T3-t)/T1)^2*(-1/T1)*(-1/T1)+6*5*4*((T3-t)/T1)^3*(-1/T1)*(-1/T1));
elseif (t>T3)&&(t<=T4)
alpha=0;
Dalpha=0.0;
D2alpha=0.0;
end
Best wishes
Torsten.
采纳的回答
Abraham Boayue
2018-5-17
编辑:Abraham Boayue
2018-5-17
Thanks for clearly formulating your problem;this is what drew my attention. It is very difficult to find a well stated problem as such. I hope my solution will be of some help.
Here is a code that I wrote to solve this problem.The code is written as if you were to solve the problem by hand.The procedure is to first implement the piecewise function and find its derivative by hand. Next, reduce each one of the second order equations into two first order equations. The end result will be a system of 4 1st order equations.For instance, I set z1 = beta and z2 = beta's derivative so that the derivative of z1 = z2 and the derivative of z2 = the double derivative of beta.This is the standard method of reducing 2nd order ode into 1st order ode. See the code below. First test the piecewise function part. It should work fine.
N = 200;
T1 = 2;
T2 = 4;
T3 = 6;
T4 = 7;
alpha_max = 30;
am = alpha_max;
t = 0 :T4/(N-1):T4;
% Implement the piecewise function
y1 = (10/T1^3)*t.^3-(15/T1^4)*t.^4+(6/T1^5)*t.^5;
y1 = am*y1;
y2 = am;
y3 = (10/T1^3)*(T3-t).^3-(15/T1^4)*(T3-t).^4+(6/T1^5)*(T3-t).^5;
y3 = am*y3;
y4 = 0;
alpha = y1.*(0<t & t<T1) +y2.*(T1<t & t<T2)+y3.*(T2<t & t<T3)+...
y4.*(T1<t & t<T2);
% Implement the derivative
y5 = (30/T1^3)*t.^2-(60/T1^4)*t.^3+(30/T1^5)*t.^4;
y5 = am*y5;
y6 = 0;
y7 = (-30/T1^3)*(T3-t).^2+(60/T1^4)*(T3-t).^3-(30/T1^5)*(T3-t).^4;
y7 = am*y7;
y8 = 0;
alpha_prime = y5.*(0<t & t<T1) +y6.*(T1<t & t<T2)+y7.*(T2<t & t<T3)+...
y8.*(T1<t & t<T2);
figure
plot(t ,alpha,'linewidth',2)
hold on
plot(t ,alpha_prime,'linewidth',2,'color','r')
legend('\alpha','\alpha_{prime}');
a = title('\alpha(t)');
set(a,'fontsize',14);
a = ylabel('y');
set(a,'Fontsize',14);
a = xlabel('t [0 7]');
set(a,'Fontsize',14);
xlim([0 T4])
grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Now run the main code:
function Twosecond_oder_odes
% Verify that your equations are implemented correctly.
% Run these functions as a single m-file.
N = 200;
T = 12;
t = 0 :T/(N-1):T;
ic = [.5 .5 1 0]; % Initial conditions Must be 4 numbers, try different values
[t,z]= ode45(@RHS, t,ic );
figure
plot(t ,z(:,1),'linewidth',2,'color','b')
hold on
plot(t ,z(:,2),'linewidth',2,'color','r')
plot(t ,z(:,3),'linewidth',2,'color','m')
plot(t ,z(:,4),'linewidth',2,'color','g')
legend('\beta','\beta_{prime}','x',',x_{prime}');
a = title('\beta(t) x(t)');
set(a,'fontsize',14);
a = ylabel('y');
set(a,'Fontsize',14);
a = xlabel('t [0 7]');
set(a,'Fontsize',14);
Ylim([-18000 16000])
grid
function dZdt= RHS(t,z)
% parameters
J2 = 0.90;
J3 = 0.10;
mA = 10;
mB = 35;
r1 = 0.25;
r2 = 0.55;
r3 = 0.15;
cc = 300;
c = 400;
kc = 2000;
k = 6000;
T1 = 2;
T2 = 4;
T3 = 6;
alpha_max = 30;
am = alpha_max;
% Coefficients from Equation 1
a1 = J2 + 2*J3 + mA*r3^2;
a2 = 2*cc*r2^2 + c*r3^2;
a3 = -c*r3;
a4 = 2*kc*r2^2 + k*r3^2;
a5 = -k*r3;
% Coefficients from Equation 2
b1 = mB;
b2 = -c*r3;
b3 = c;
b4 = a5;
b5 = k;
% Implements the piecewise function, alpha (Test this to make sure its
% right)
y1 = (10/T1^3)*t.^3-(15/T1^4)*t.^4+(6/T1^5)*t.^5;
y1 = am*y1;
y2 = am;
y3 = (10/T1^3)*(T3-t).^3-(15/T1^4)*(T3-t).^4+(6/T1^5)*(T3-t).^5;
y3 = am*y3;
y4 = 0;
alpha = y1.*(0<t & t<T1) +y2.*(T1<t & t<T2)+y3.*(T2<t & t<T3)+...
y4.*(T1<t & t<T2);
% Implements the derivative of alpha (Do it manually!)
y5 = (30/T1^3)*t.^2-(60/T1^4)*t.^3+(30/T1^5)*t.^4;
y5 = am*y5;
y6 = 0;
y7 = (-30/T1^3)*(T3-t).^2+(60/T1^4)*(T3-t).^3-(30/T1^5)*(T3-t).^4;
y7 = am*y7;
y8 = 0;
alpha_prime = y5.*(0<t & t<T1) +y6.*(T1<t & t<T2)+y7.*(T2<t & t<T3)+...
y8.*(T1<t & t<T2);
% This is the right hand side of Equation 1
F = 2*cc*r1*r2*alpha_prime + 2*kc*r1*r2*alpha;
% Implements two 2nd orders into 4 1st oderes ODES
dZdt_1 = z(1); % Solution for beta (z1 = beta, z2 = derivative of beta)
dZdt_3 = z(3); % Solution for x (z3 = x, z4 = derivative of x)
dZdt_2 = a2*z(2)+a4*z(1)+a3*z(4)+a5*z(3)-F; % dZdt_2 = double derivative of beta
dZdt_2 = (-1/a1)*dZdt_2;
dZdt_4 = b2*z(2)+b4*z(1)+b3*z(4)+b5*z(3);% dZdt_4 = double derivative of x
dZdt_4 = (-1/b1)*dZdt_4;
dZdt =[dZdt_1; dZdt_2; dZdt_3;dZdt_4];
更多回答(1 个)
Walter Roberson
2018-5-17
None of the ode* routines can handle piecewise expressions. All of them require that the calculation in any one ode* call is continuous and that the first two derivatives of everything you programmed is also continuous (numeric estimates of them are created.) If the first derivative is not continuous then often a message about integration tolerance is shownń if it is the second derivatives that are not then it will typically just produce an incorrect solution.
You need to recode so that every time that your piecewise would switch branches, that you fire an event function set to terminal. Then call ode* again for the remainder of the tspan, using the last output of the previous run as the initial condition of the next. Keep doing this until the full tspan is processed.
There was a case about 2 weeks ago where the event function was firing every 3 to 7 calls. I posted example code there showing how to run the ode.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!