ODE45 not updating a variable of nonlinear 2nd order diff.eq system

3 ビュー (過去 30 日間)
Res
Res 2018 年 9 月 8 日
回答済み: Res 2018 年 9 月 11 日
Hello everyone, I was trying to solve a simple nonlinear system of diff.eqn (piece-wise linear) that involves coulomb friction. The problem I face is x(3) does not update after passing if condition... Although when i debug it, I can see that its varying accordingly... but at the end get a vector filled with a constant number(initial condition x3(0).) Could you please guide me how to cope with this issue???
function [dxdt]=func6(t,x)
C=10;
k=5;
kd=10;
F=25;
mu=0.5;
Nf=15;
m=1;
wn=sqrt(k/m);
w=wn;
fE=F*sin(w*t);
dxdt=zeros(size(x));
x1=x(1);
x2=x(2);
x3=x(3);
dxdt(1)=x(2);
%
nonlinear force calculation function
F_NL=kd*(x(1)-x(3));
if F_NL>=mu*Nf
F_NL=mu*Nf;
x(3)=x(1)-(mu*Nf)/kd;
end
if F_NL<=-mu*Nf
F_NL=-mu*Nf;
x(3)=x(1)+(mu*Nf)/kd;
end
dxdt(2)=1/m*(fE-F_NL-C*x(2)-k*x(1));
end
%%%and am using this script to run my ode calculations...
h=0.01;
tspan = 0:h:10;
initialvalues=[0 0 0];
[t,x]=ode23(@func6,tspan,initialvalues);

採用された回答

Aquatris
Aquatris 2018 年 9 月 9 日
編集済み: Aquatris 2018 年 9 月 9 日
The reason is simple. Your dxdt(3) is always equal to 0. You do not have an equation for dxdt(3) in your function other than the first zeros(size(x)) assignment. Since the initial value is also zero, x(3) stays 0 at all times.
By any chance, those if loops supposed to be for dxdt(3) not x(3)?
  4 件のコメント
Aquatris
Aquatris 2018 年 9 月 11 日
Something like this, it needs some modifications though;
The function;
function [dxdt]=func6(t,x)
% specify the global variables
% the function needs to use
global x3 counter t_x3
C=10;
k=5;
kd=10;
F=25;
mu=0.5;
Nf=15;
m=1;
wn=sqrt(k/m);
w=wn;
fE=F*sin(w*t);
dxdt=zeros(size(x));
x1=x(1);
x2=x(2);
dxdt(1)=x(2);
nonlinear force calculation function
F_NL=kd*(x(1)-x3(counter));
counter = counter +1;% counter to store x3 values
t_x3(counter) = t; % time for variable x3, ode is weird sometimes
% and the size of the time output of ODE did
% not match the x3 variable
if F_NL>=mu*Nf
F_NL=mu*Nf;
x3(counter)=x(1)-(mu*Nf)/kd;
elseif F_NL<=-mu*Nf
F_NL=-mu*Nf;
x3(counter)= x(1)+(mu*Nf)/kd;
else
x3(counter) = 0; % I assigned 0 but what happens here
% is up to your problem
end
dxdt(2)=1/m*(fE-F_NL-C*x(2)-k*x(1));
end
And the script;
clear,clc
% specify the global variables
global x3 counter t_x3
x3 = 0;
counter = 1;
h=0.01;
tspan = 0:h:10;
initialvalues=[0 0 ]';
[t,x]=ode45(@func6,tspan,initialvalues);
plot(t,x,t_x3,x3,'.')

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

その他の回答 (1 件)

Res
Res 2018 年 9 月 11 日
Thank you, really appreciated.

Community Treasure Hunt

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

Start Hunting!

Translated by