抛物型方程求解

相关主题

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

10,01,01(,0),01(0,),(1,),01

(,)x t t x t

u u x t t x u x e x u t e u t e t u x t e ++∂∂-=<<<≤∂∂=≤≤==<≤=

运行：前向euler 法

[xx,tt,uh]=equationepaowu2('myfun','myfun1','myfun1','myfun2',1,[0,1],[0,1],[1/10,1/200]); function [xx,tt,uh]=equationepaowu2(myfun,myfun1,myfun2,myfun3,a,xxx,ttt,step) %利用差分方法求抛物型方程数值解；

%myfun--方程右端f(x,t);

%myfun1--u(x,0);

%myfun2--u(t1,t);

%myfun3--u(t2,t);

%[x1,x2]--x 的取值范围；

%[t1,t2]--t 的取值范围；

%a-正常数

%h,tao-分别是x,t 方向的步长。

%——————————————————————

%激活函数

f=fcnchk(myfun);

f1=fcnchk(myfun1);

f2=fcnchk(myfun2);

f3=fcnchk(myfun3);

x1=xxx(1);x2=xxx(2);

t1=ttt(1);t2=ttt(2);

h=step(1);tao=step(2);

%__________________________________

%划分网格,x1-nt+1行，nx+1列。

x=linspace(x1,x2,round((x2-x1)/h)+1);

t=linspace(t1,t2,round((t2-t1)/tao)+1);

nx=size(x,2);

nt=size(t,2);

[xx,tt]=meshgrid(x,t);

%________________________________________

%赋初值及边值

size(x1)

size(x)

U0=zeros(size(xx));

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

U0(2:nt,1)=f2(t(2:nt))';

U0(2:nt,nx)=f3(t(2:nt))';

F=f(xx,tt);

%_________________________________

%计算步长

r=(a*tao)/h^2

%_______________________

%迭代；由低层向高层计算,i-代表层数；

for i=2:nt

for j=2:nx-1

ss=(1-2*r)*U0(i-1,j)+r*(U0(i-1,j-1)+U0(i-1,j+1))+tao*F(i-1,j);

U0(i,j)=ss;

end

uh=U0;

%作图

uh=U0;

Uh=exp(xx+tt);

ee=max(max(abs(Uh-uh)))

mesh(xx,tt,Uh-uh);

function F=myfun(x,t)

F=0*x.*t;

function f1=myfun1(x)

f1=exp(x)

function f1=myfun2(x)

f1=exp(1+x)

后向euler法

function [xx,tt,uh]=equationepaowu222(myfun,myfun1,myfun2,myfun3,a,xxx,ttt,step) %利用差分方法求抛物型方程数值解；

%myfun--方程右端f(x,t);

%myfun1--u(x,0);

%myfun2--u(t1,t);

%myfun3--u(t2,t);

%[x1,x2]--x的取值范围；

%[t1,t2]--t的取值范围；

%a-正常数

%h,tao-分别是x,t方向的步长。

%——————————————————————%激活函数