一维波动方程数值解 Matlab 教程（从入门到出图）

#一维波动方程数值解 Matlab 教程（从入门到出图）| 来源: 网络整理| 查看: 265

最近在计算波动方程数值解时，发现网上很难找到比较详细实用的Matlab教程，所以我打算自己写一个，以一维弦的波动为例，从理论推导到代码实现，跟大家分享一下。

杜帅：一维波动方程数值解 Matlab 教程（从入门到出图）——1波动方程基本概念杜帅：一维波动方程数值解 Matlab 教程（从入门到出图）——2差分方程基本概念杜帅：一维波动方程数值解 Matlab 教程（从入门到出图）——3数值计算的Matlab实现

前两篇已经讲了一维波动方程以及差分方法数值计算的基本知识，最后我们来看看如何用Matlab实现二阶偏微分方程的数值计算[1]：

一、一维波动方程的差分形式

微分方程

$c^2\frac{\partial^2 u}{\partial x^2}= \frac{\partial^2 u}{\partial t^2}\\$

我们将自变量 x 和 t 离散化，分别均分为 N_x 和 N_t 份，步长分别为 h_x 和 h_t ，其中第 i 个 x 和第 j 个 t 对应的位移为 $u_{i,j}=u(x_i,t_j)$

根据二阶差商的写法，改写为差分形式

$c^2\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h_x^2}=\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_t^2}\\$

整理得

$u_{i,j+1}=(c\frac{h_t}{h_x})^2(u_{i+1,j}-2u_{i,j}+u_{i-1,j})+2u_{i,j}-u_{i,j-1}\\$

这样等号右边都是关于时间 t_j 和 $t_{j-1}$ 的位移，只有左边是关于 $t_{j+1}$ 的。

而对于位置 $x_{i-1}$ 、 x_i 和 $x_{i+1}$ ，初始条件里已经全部给出，所以只需要带入上式就可以循环计算出所有结果。

边界条件的差分形式

位移边界条件 u(0,t)=u(l,t)=0

即 $u_{1,j}=0\\ u_{N_x,j}=0$

应力边界条件 $EA\frac{\partial u(l,t)}{\partial x}=F$

差分形式 $\frac{u_{N_x,j}-u_{N_x -1,j}}{h_x}=\frac{F}{EA}$

即 $u_{N_x,j}=u_{N_x -1,j}+\frac{Fh_x}{EA}\\$

二、Matlab程序实现

1、参数离散化

c = 1; % 波速 xmin = 0; xmax = 10; Nx = 1000; % 空间离散化 tmin = 0; tmax = 15; Nt = 2000; % 时间离散化 x = linspace(xmin, xmax, Nx)'; t = linspace(tmin, tmax, Nt); hx = (xmax - xmin)/(Nx - 1); ht = (tmax - tmin)/(Nt - 1);

2、初始条件

我们设置初始位移和初始速度均为零。

u = zeros(Nx, Nt);

3、循环计算

只需要从 j=2 开始循环时间，带入初始条件，计算下一时刻所有位置的位移。

其中二次差分可以利用Matlab的向量运算，避免再套循环。由于二次差分后少了两个元素，我们在两端补零即可，因为后边边界条件还会对两端进行修正，所以没有影响。

for j = 2:Nt-1 un = u(:,j); ddu = [0; un(1:end-2) - 2*un(2:end-1) + un(3:end); 0]; u(:,j+1) =(c*ht/hx)^2 * ddu+ 2*u(:,j) - u(:,j-1); u(1, j+1) = 0;%左端位移边界条件 %右端应力边界条件 if t(j)

【本文地址】

公司简介

联系我们