计算科学 - 使用有限差分法MATLAB的一维波动方程 - 吾爱随笔录

使用有限差分法MATLAB的一维波动方程

计算科学 pde 有限差分边界条件

2021-12-06 14:00:47

我有波动方程

u_{t t} = 4 u_{x x}

$u_{tt} = 4 u_{xx}$ 有边界条件

u (0, t) = u (L, t) = 0, x \leq 0 \leq 2 π, t \geq 0

$u(0,t) = u(L,t) = 0\,,\quad x \leq 0 \leq 2\pi \,,\quad t\geq 0$ 和初始条件

\begin{aligned} u (x, 0) = {\begin{cases} 1 & for π - 1 \leq x \leq π + 1 \\ 0 & otherwise \end{cases} \\ u_{t} (x, 0) = x^{2} s i n^{2} (4 x) \end{aligned}

$\begin{align} &u(x,0)=\begin{cases}1 & \text{for} \quad \pi-1 \leq x \leq \pi +1 \\0 &\text{otherwise} \end{cases} \\[0.5em] &u_{t}(x,0) = x^2sin^2(4x) \end{align}$

我试图通过有限差分法在 MATLAB 上实现这个问题，并使用该surf函数将其绘制为 3D 波；但是，我遇到的问题是如何编写第一个初始条件。

我正在使用的功能是

function [u, q] = Wave(f1,f2,g0,g1,xspan,tspan,nx,nt,alpha)
x0 = xspan(1)
xf = xspan(2)
t0 = tspan(1)
tf = tspan(2)
dx = (xf - x0)/nx;
dt = (tf-t0)/nt;
x = [0:nx]'*dx;
t = [0:nt]*dt;
q = alpha*(dt/dx)^2;
q1 = q/2;
q2 = 2*(1-q);
u(:,1) = f1(x);
for k = 1:nt+1, u([1 nx+1],k) = [g0(t(k)); g1(t(k))]; end
u(2:nx,2) = q1*u(1:nx-1,1) + (1-q)*u(2:nx,1) + q1*u(3:nx+1,1) + dt*f2(x(2:nx));
for k = 3:nt+1
    u(2:nx,k) = q*u(1:nx-1,k-1) + q2*u(2:nx,k-1) + q*u(3:nx+1,k-1) - u(2:nx,k-2);
end

surf(t,x,u) 
xlabel('t')
ylabel('x')
zlabel('u(x,t)')
end

我知道在函数中，变量f1是由第一个初始条件控制的。在函数中，我不知道如何将初始条件合并到方法中。我假设要使用 for 循环和 if 语句，但我尝试的任何方法都不起作用。

我用来绘制图表的脚本是

clc
clear all
f1 = @(x) ////;
f2 = @(x) (x.^2).*((sin(4.*x)).^2);
g0 = @(t) 0;
g1 = @(t) 0;
xspan = [0,2*pi];
tspan = [0,1];
nx = 20;
nt = 40;
alpha = 4;
[u,r] = Wave(f1,f2,g0,g1,xspan,tspan,nx,nt,alpha);

正如在脚本中看到的那样f1 = @(x) ////;。我不确定在初始条件下放什么。任何帮助将不胜感激。

1个回答

检查 x 的间隔，然后将逻辑转换为 double。

f1 = @(x) double(pi-1 <= x && x <= pi+1);

或者如果初始条件不是1和0，写一个函数

function u = f1(x)
    u = zeros(size(x));
    is_one = pi-1 <= x && x <= pi+1;
    u(is_one) = 1;
    % u(~is_one) = 0;  not required, just for your reference
end

编辑：

$u(x, 0)$ 由

u(1,:) = f1(x);

然后你用它来计算 $u(x,\Delta t)$ ，即u(:,2)在这一行

u(2:nx,2) = q1*u(1:nx-1,1) + (1-q)*u(2:nx,1) + q1*u(3:nx+1,1)...
            + dt*f2(x(2:nx));

你已经做到了，不是吗？

其它你可能感兴趣的问题

上一篇编译不同的 C 源文件时生成相同的 .so 共享库（没有符号）下一篇如何绘制彩虹？