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

I=10;???%x方向分段數(shù)
dx=0.2;?%x方向的步長
T=0.2;??%計算總時間
N=c;????%t方向的分段數(shù)
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））??%當(dāng)精度不滿足時，繼續(xù)迭代
????h0=h;??%保存前一次計算的h計算結(jié)果
????v0=v;??%保存前一次計算的v計算結(jié)果
????for?n=2:N+1???%計算h和v需要用內(nèi)部點表示的邊界
????????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

????%計算前后兩個矩陣對應(yīng)元素的最大誤差
????eh=max（max（abs（（h0-h）./（h+eps））））;???
????ev=max（max（abs（（v0-v）./（v+eps））））;
????NN=NN+1;????%記錄循環(huán)次數(shù)
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方向的網(wǎng)格線，以便繪制三維圖
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實現(xiàn)偏微分方程的差分計算?源程序代碼\chafenfangcheng.m

?????文件????????294??2014-02-06?11:28??MATLAB實現(xiàn)偏微分方程的差分計算?源程序代碼\聯(lián)系我們.txt

?????目錄??????????0??2014-02-06?11:37??MATLAB實現(xiàn)偏微分方程的差分計算?源程序代碼

-----------?---------??----------?-----??----

?????????????????2192????????????????????3