1

我正在研究一个项目，该项目与使用具有以下边界条件的 MATLAB 中的中心差近似以数值方式求解 2D (x, y, t) 波动方程有关：

阴谋

总装公式为：

FD公式

我了解一些边界条件（BC），例如

在 j=m 处，du/dy=0，

公元前,

但我不确定如何在 MATLAB 中实现这些边界条件。

一个朋友给了我这些方程式：

情商

这是我对 MATLAB 代码的尝试，
但我无法进一步取得进展：

% The wave function
% Explicit 

% Universal boundary conditions for all 3 cases:
% u=0 at t=0
% du/dt=0 at t=0

% Case 1 boundary conditions
% At x=0, u=2sin(2*pi*t/5);
% At y=0, du/dy=0;
% At y=2, du/dy=0;
% At x=5, du/dx=0;
% u=0 and du/dt=0 at t=0;
%-------------------------------------------------------------------------%


% Setting up
clc; clear all; close all;

% length, time, height
L = 5; % [m]
h = 2; % [m]
T = 10; % [s]

% Constants
c_x = 1; % arbitrary 
c_y = 1; % arbitrary

dx = 0.1; % <x> increment
dy = 0.1; % <y> increment
dt = 0.1; % time increment
nx = L/dx + 1; % number of <x> samples
ny = h/dy + 1; % number of <y> samples
nt = T/dt + 1; % number of time samples
t(:,1) = linspace(0, T, nt); 

theta_x = c_x*(dt^2)/(dx^2); 
theta_y = c_y*(dt^2)/(dy^2); 
% theta_x = theta_y
theta = theta_x;
%-------------------------------------------------------------------------%

% The matrix 
U = zeros(nt, nx, ny);

% Setting up the <U> matrix with the boundary conditions - case 1
U(1, :, :) = 0; % U=0 at t=0

for tt=1:nt   % U=2sin(2pi/5*t) at x=0
    for jj=1:ny
        U(tt, 1, jj)=2*sin(2*pi/5.*t(tt));
    end 
end


for it=2:t
    for ix=2:nx-1
        for iy=2:ny-1
            % Boundary conditions

            % General case (internal):
            U1 = -U(it-1, ix, iy);
            U2 = 2*(1-2*theta)*u(it, ix, iy);
            U3 = theta*U(it, ix-1, iy);
            U4 = theta*U(it, ix+1, iy);
            U5 = theta*U(it, ix, iy-1);
            U6 = theta*U(it, ix, iy+1);



        end 
    end 
end
4

1 回答 1

0

您拥有的总装公式也适用于边界。j = 1复杂之处在于，当您在和时应用公式时j = m，您的j = 0和j = m+1术语不在您的网格中。

为了改善这个问题，边界条件为您提供了网格外点和网格上点之间的关系。

正如你所指出的，dudy = 0条件给了你u(i,m-1) == u(u,m+1)边界上的关系。因此，您使用通用汇编公式并将所有m+1术语替换m-1 为 边界。对于下边界，您也会有类似的关系。

于 2015-03-12T21:01:49.713 回答