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

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.

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}=f\left(x,t\right)$ hay viết ở dạng khác ${\dot{x}}\left(t\right)=f\left(x,t\right)$
Ví dụ: để giải phương trình vi phân bậc nhất ${\dot{x}}\left(t\right)=3e^{-t}$ với điều kiện ban đầu $x\left(0\right)=0$ bằng Matlab Ode45

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(t\right)+5\dot{x}\left(t\right)-4x\left(t\right)=sin\left(10t\right)
\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}=x\left(t\right)\\
x_{2}=\dot{x}\left(t\right)
\end{matrix}\right.
\Rightarrow
\left\{\begin{matrix}
\dot{x}_{1}=\dot{x}\left(t\right)\\
\dot{x}_{2}=\ddot{x}\left(t\right)
\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}+sin\left(10t\right)
\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\\
sin\left ( 10t \right )
\end{bmatrix}
\end{equation}

Giải bằng Matlab

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}
M\ddot{x}\left(t\right)+c\dot{x}\left(t\right)+kx\left(t\right)=F\left(t\right)
\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}=x\left(t\right)\\
x_{2}=\dot{x}\left(t\right)
\end{matrix}\right.
\Rightarrow
\left\{\begin{matrix}
\dot{x}_{1}=\dot{x}\left(t\right)\\
\dot{x}_{2}=\ddot{x}\left(t\right)
\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{F\left(t\right)}{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}F\left(t\right)
\end{equation}
Code Matlab

Kết quả
ex4

Tiến Anh

Xin chào, tôi là Tiến Anh. Nội dung trên blog này là những chia sẽ kiến thức của tôi với hy vọng sẽ hữu ích cho mọi người.

3 Responses

  1. NTP PRO viết:

    Cảm ơn bài viết rất hay của bạn! Mình có 1 góp ý nhỏ ở phần “Giải ptvp bậc 2” là: hệ phương trình nhận được phải là dx1/dt = x2; dx2/dt = -5*x2 + 4*x1 + sin(10*t).

  2. Đỗ Tùng Lâm viết:

    cảm ơn bạn nhé!

  3. Mai Anh viết:

    cảm on bạn rất nhiều

Trả lời