有限元法解偏微分方程程序

有限元法解偏微分方程程序

clear;clc;close all
c=1000;
g=9.81;
lambda=0.025;
delta=0.005;

I=10;   %x方向分段数
dx=0.2; %x方向的步长
T=0.2;  %计算总时间
N=c;    %t方向的分段数
dt=dx/c; %t方向的步长

A=1000;
w=1000;
%定义计算矩阵h和v
h=zeros（N+1I+1）;
v=zeros（N+1I+1）;

%计算h和v能确定的两个边界
%初始条件
h（1:）=A;
v（1:）=0;
t=0:dt:T;
%边界条件
h（:1）=A*cos（w*t）‘;
v（:I+1）=0;
NN=0;
eh=1;  %定义误差
ev=1;  %定义误差

msgbox（‘MATLAB编程答疑，请加QQ: 1530497909‘‘MATLAB答疑‘‘help‘）
web http://url.cn/TKcdXk -browser

while（（ev>1e-5）|（eh>1e-5））  %当精度不满足时，继续迭代
    h0=h;  %保存前一次计算的h计算结果
    v0=v;  %保存前一次计算的v计算结果
    for n=2:N+1   %计算h和v需要用内部点表示的边界
        v（n1）=（g/c）*（h（n1）-h（n-12））+v（n-12）*（1-abs（v（n-12））*lambda*dt/（4*delta））;
        h（nI+1）=h（n-1I）+（c/g）*v（n-1I）-v（n-1I）*abs（v（n-1I））*lambda*dx/（4*g*delta）;
    end
    %计算中间点 
    for i=2:I
        for n=2:N+1
            h（ni）=0.5*（（h（n-1i-1）+h（n-1i+1））+（c/g）*（v（n-1i-1）-v（n-1i+1））-（lambda*dx/（4*g*delta））*（v（n-1i-1）*abs（v（n-1i-1））-v（n-1i+1）*abs（v（n-1i+1））））;
            v（ni）=0.5*（（g/c）*（h（n-1i-1）-h（n-1i+1））+（v（n-1i-1）+v（n-1i+1））-（lambda*dt/（4*delta））*（v（n-1i-1）*abs（v（n-1i-1））+v（n-1i+1）*abs（v（n-1i+1））））;
        end
    end

    %计算前后两个矩阵对应元素的最大误差
    eh=max（max（abs（（h0-h）./（h+eps））））;   
    ev=max（max（abs（（v0-v）./（v+eps））））;
    NN=NN+1;    %记录循环次数
end

figure
subplot（311）
plot（th（:1））
title（‘i=0处的h值随时间变化曲线‘）

subplot（312）
plot（th（:6））
title（‘i=5处的h值随时间变化曲线‘）

subplot（313）
plot（th（:11））
title（‘i=10处的h值随时间变化曲线‘）

[xxtt]=meshgrid（0:It）;   %绘制x和y方向的网格线，以便绘制三维图
figure
mesh（xxtth）    %绘制h
xlabel（‘x‘）
ylabel（‘t‘）
zlabel（‘h‘）
figure
mesh（xxttv）   %绘制v
xlabel（‘x‘）
ylabel（‘t‘）
zlabel（‘v‘）

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       1898  2014-02-06 11:21  MATLAB实现偏微分方程的差分计算 源程序代码\chafenfangcheng.m

     文件        294  2014-02-06 11:28  MATLAB实现偏微分方程的差分计算 源程序代码\联系我们.txt

     目录          0  2014-02-06 11:37  MATLAB实现偏微分方程的差分计算 源程序代码

----------- ---------  ---------- -----  ----

                 2192                    3