Kiến thức

Giải phương trình vi phân bằng Matlab ode45

Giải phương trình vi phân bằng Matlab ode45

September 11, 2016

0

Bài hướng dẫn này sẽ giới thiệu các bạn cách giải phương trình vi phân bậc 2 có hệ số là các hằng số bằng Matlab ode45. Ode là một công cụ được mình sử dụng khá thường xuyên để giải các phương trình vi phân trong matlab. Chúng ta cần nhớ rằng ode45 chỉ có thể giải được phương trình vi phân bậc nhất. Vì vậy các bạn muốn giải phương trình vi phân bậc thì chúng ta phải chuyển về hệ hai phương trình vi phân bậc nhất.

Bạn đang xem: Giải phương trình vi phân bằng Matlab ode45

Phương trình vi phân bậc nhất

Phương trình vi phân bậc nhất có dạng: $frac{dx}{dt}=fleft(x,tright)$ hay viết ở dạng khác ${dot{x}}left(tright)=fleft(x,tright)$
Ví dụ: để giải phương trình vi phân bậc nhất ${dot{x}}left(tright)=3e^{-t}$ với điều kiện ban đầu $xleft(0right)=0$ bằng Matlab Ode45

function first_oder_ode % SOLVE dx/dt = -3 exp(-t). % initial conditions: x(0) = 0 t=0:0.001:5; % time scalex initial_x=0; [t,x]=ode45( @rhs, t, initial_x); plot(t,x); xlabel('t'); ylabel('x'); grid on function dxdt=rhs(t,x) dxdt = 3*exp(-t); end end 

Kết quả:
ex1

Giải phương trình vi phân bậc hai

Ví dụ: chúng ta tiến hành giải phương trình vi phân bậc hai sau
begin{equation}label{ref1}
ddot{x}left(tright)+5dot{x}left(tright)-4xleft(tright)=sinleft(10tright)
end{equation}
Như đã nói ở trên, chúng ta chuyển phương trình vi phân bậc 2 thành hệ hai phương trình vi phân bậc nhất bằng cách đặt ẩn phụ $x_{1}, x_{2}$:
$left{begin{matrix}
x_{1}=xleft(tright)\
x_{2}=dot{x}left(tright)
end{matrix}right.
Rightarrow
left{begin{matrix}
dot{x}_{1}=dot{x}left(tright)\
dot{x}_{2}=ddot{x}left(tright)
end{matrix}right.$
Thế vào phương trình eqref{ref1} chúng ta có hệ
begin{aligned}
left{begin{matrix}
x_{1}=&x_{2}\
x_{2}=&-5x_{2}+4x_{1}+sinleft(10tright)
end{matrix}right.
end{aligned}
Nếu đặt $X=begin{bmatrix}
x_1\
x_2
end{bmatrix}$ thì hệ trên có thể viết thành dạng phương trình vi phân bậc nhất của biến $X$
begin{equation}
X=begin{bmatrix}
0 & 1\
4 & -5
end{bmatrix}X+begin{bmatrix}
0\
sinleft ( 10t right )
end{bmatrix}
end{equation}

Giải bằng Matlab

function second_oder_ode % SOLVE d2x/dt2+5 dx/dt - 4 x = sin(10 t) % initial conditions: x(0) = 0, x'(0)=0 t=0:0.001:3; % time scale initial_x = 0; initial_dxdt = 0; [t,x]=ode45( @rhs, t, [initial_x initial_dxdt] ); plot(t,x(:,1)); xlabel('t'); ylabel('x'); grid on function dxdt=rhs(t,x) dxdt_1 = x(2); dxdt_2 = -5*x(2) + 4*x(1) + sin(10*t); dxdt=[dxdt_1; dxdt_2]; end end

Kết quả
ex2

Ứng dụng trong mô phỏng

Dưới đây sẽ là ví dụ dùng ode45 giải phương trình vi phân trong mô phỏng hệ vật nặng, lò xo và giảm chấn
ex3
Hệ lò xo giảm chấn có 1 bậc tự do, hệ thống bao gồm xe có khối lượng $M$, gắn với lò xo độ cứng $k$ và giảm chấn có hệ số $c$. Khi ngoại lực $F(t)$ tác động lên xe thì chúng ta cần xác định được vị trí $x(t)$ của xe. Áp dụng định luật 2 Newton chúng ta xác định được phương trình động lực học của cơ hệ
begin{equation}label{ref2}
Mddot{x}left(tright)+cdot{x}left(tright)+kxleft(tright)=Fleft(tright)
end{equation}
Viết lại ở dạng hệ phương trình vi phân bậc nhất bằng cách đặt ẩn phụ $x_{1}, x_{2}$:
begin{equation*}left{begin{matrix}
x_{1}=xleft(tright)\
x_{2}=dot{x}left(tright)
end{matrix}right.
Rightarrow
left{begin{matrix}
dot{x}_{1}=dot{x}left(tright)\
dot{x}_{2}=ddot{x}left(tright)
end{matrix}right.end{equation*}
Thế vào eqref{ref2} chúng ta được hệ phương trình vi phân bậc nhất
begin{aligned}
left{begin{matrix}
dot{x}_1=&x_2\
dot{x}_2=&-dfrac{c}{M}x_2-dfrac{k}{M}x_1+dfrac{Fleft(tright)}{M}
end{matrix}right.
end{aligned}
Hay dạng phương trình trạng thái của cơ hệ
begin{equation}
begin{bmatrix}
dot{x}_1\
dot{x}_2
end{bmatrix}=begin{bmatrix}
0 & 1\
-dfrac{k}{M} & -dfrac{c}{M}
end{bmatrix}begin{bmatrix}
x_1\
x_2
end{bmatrix}+begin{bmatrix}
0\
dfrac{1}{M}
end{bmatrix}Fleft(tright)
end{equation}
Code Matlab

function mass_spring_damper() clear; close all; t_start = 0; t_end = 6; %final time in seconds. time_span =t_start:0.001:t_end; k = 40; % spring stiffness. N/m m = 5; % mass, kg c = 31; initial_position = 0; initial_speed = 0; x0 = [initial_position initial_speed]; [t,x]=ode45(@rhs,time_span,x0); plot(t,x(:,1)); ylim([-.1 .5]); grid on %************************************** % solves m x''+ c x' + k x = f(t) %************************************** function xdot=rhs(t,x) xdot_1 = x(2); xdot_2 = -(c/m)*x(2) - (k/m)*x(1) + force(t)/m; xdot = [xdot_1 ; xdot_2 ]; end %******************** % The forcing function, edit to change as needed. %******************** function f=force(t) P = 100; % force amplitude %f=P*sin(omega*t); f=10; %unit step %if t<eps %impulse % f=1 %else % f=0; %end %f=P*t; %ramp input end end

Kết quả
ex4

Tags:

Matlab ode45

Chuyên mục: Kiến thức

Related Articles

Trả lời

Email của bạn sẽ không được hiển thị công khai. Các trường bắt buộc được đánh dấu *

Check Also
Close
Back to top button