偏微分方程报告

2009级数学与应用数学和信息与计算科学专业

偏微分方程数值解上机实验

实验题目 利用有限元方法和有限差分方法求解偏微分方程 完成日期 2012年12月17日 学生姓名 张灵刚 所在班级 1102090 任课教师 王晓东

西北工业大学理学院应用数学系

目录

一．实验目的………………………………………..…….（2） 二．实验要求…………………………………….…...….（2） 三．实验题目……………………………………….………（3） 四．实验二……………………………………………….…（4）

1.实验内容…………………………………….…..（4） 2.实验原理……………………………………….…（4）. 3算法流程………………………………….………（5） 4结果分析…………………………………………..（5）

5总结讨论…………………………………………..（6） 6源程序………………………………………..……（6）

五. 实验三

1.实验内容…………………………………………（17） 2.实验原理…………………………………….….（17）. 3算法流程………………………………………….（18） 4结果分析……………………………………….….（18）.

5总结讨论………………………………………….（21） 6源程序……………………………………………（21）

偏微分方程数值解上机实验报告 实验地点：数学系机房

实验时间：第13—15周，周一、四下午5、6节 实验分数：占期末考试成绩的30% 一、实验目的及意义

掌握有限元方法和有限差分方法的程序实现；学会选择合适的有限差分格式求解一维非线性对流占优的非定常对流扩散问题；学会使用三角线性元和四边形线性元的有限元方法求解二维椭圆方程边值问题，并对计算结果进行收敛性分析；尝试采用有限元方法或有限差分方法实现二维初边值抛物型方程的大规模数值求解。通过实验可以提高学生的动手能力，加深学生对算法的理解。 二、实验要求

在下列给出的三个问题中，最少选择两个问题进行编程实现。要求给出格式的推导过程、算法流程、实现程序、选取的网格参数、以表格或图形的方式给出计算结果、对计算结果进行分析、最后对实验进行总结和讨论。

问题2：用三角线性元和四边形线性元的有限元方法求解方程：

u82cos(2x)sin(2y),(x,y)(0,1)(0,1)u(x,0)u(x,1)0;u(0,y)u(1,y)sin(2y).

取h=1/4, 1/8, 1/16, 1/32, 1/64,比较两种方法的计算精度，并给出数值收敛率.

问题3：选用合适的数值方法求解方程：

u

(44x24y2)1u,(x,y)(0,1)(0,1)t

u(0,y,t)u(1,y,t)uy(x,0,t)uy(x,1,t)0;u(x,y,0)sin(x2)cos(y2).

求t0.1,0.5,1.0时，点(

[***********]67

,)(,)(,)(,)(,)、[***********]128128

[***********]3

(,)、(,)、(,)、(,)处的数值解。 [***********]256256

上机实验（二）

一、实验内容

用三角线性元和四边形线性元的有限元方法求解方程：

u82cos(2x)sin(2y),(x,y)(0,1)(0,1)u(x,0)u(x,1)0;u(0,y)u(1,y)sin(2y).

取h=1/4, 1/8, 1/16, 1/32, 1/64,比较两种方法的计算精度，并给出数值收敛率 二、实验原理

由于计算机只能存储有限个数据和做有限次运算，所以必须把连续问题离散化，有限元法通过把求解偏微分问题转化为变分问题，实质上就是Ritz-Galerkin法，利用有穷维空间近似代替无穷维空间，从而转化为求多元二次函数的极值问题。然后选定单元形状，对求解域进行剖分。构造基函数或单元形状函数，形成有限元空间了，再形成有限元方程，并提供求解有限元方程的有效方法。 三、算法流程

四、计算结果及分析

结果分析：通过给出的上述结果可以发现三种方法相比四边形元的实际求解误差最小，这是由于四边形型元中的基函数中含有

xy项，使得四边形元的误差远比三角形元小，且当空间步长不

断加细时，向前差分的误差逐渐和四边形元接近。而对于收敛率的求解，如果舍去误差，三种方法的收敛率都约为3. 五、总结和讨论

有限元计算的有关问题有：把初值问题化为变分形式，对求解域作网络分割，构造基函数（或单元形状函数），形成有限元

方程。通过本次实验，我懂得了用有限元方法求解初值问题的基本数学思想,理解了线性有限元法的基本原理和方法。另外，我也懂得了按Galerkin方法推导有限元方程的优点，它比Ritz法更加方便直接。我也对虚功原理有了初步的认识。因为Galerkin方法基于虚功原理，所以不但可用于保守场问题，也可使用于非保守场即非驻定问题。 附录：程序源代码 1.三角形元的型函数

function [phi, dxphi, dyphi]=shape(point, element, ei, gx) x1=point(element(ei,1),1); x2=point(element(ei,2),1); x3=point(element(ei,3),1); y1=point(element(ei,1),2); y2=point(element(ei,2),2); y3=point(element(ei,3),2);

S=max(y2-y1,x2-x1)*max(y3-y1,x3-x1);

ksi=((x3*y1-y3*x1)+(y3-y1)*gx(1)+(x1-x3)*gx(2))/S;

eta=((x1*y2-y1*x2)+(y1-y2)*gx(1)+(x2-x1)*gx(2))/S; phi(1)=1-ksi-eta; phi(2)=ksi; phi(3)=eta;

dxphi(1)=-(y3-y1)/S-(y1-y2)/S; dxphi(2)=(y3-y1)/S; dxphi(3)=(y1-y2)/S;

dyphi(1)=-(x1-x3)/S-(x2-x1)/S; dyphi(2)=(x1-x3)/S; dyphi(3)=(x2-x1)/S; 2.四边形元的型函数

function [phi, dxphi, dyphi]=shape1(point, element, ei, gx)

x1=point(element(ei,1),1); x2=point(element(ei,2),1);

x3=point(element(ei,3),1); x4=point(element(ei,4),1); y1=point(element(ei,1),2); y2=point(element(ei,2),2); y3=point(element(ei,3),2); y4=point(element(ei,4),2);

S=max(y2-y1,x2-x1)*max(y4-y1,x4-x1); ksi=(gx(1)-x1)/(x2-x1); eta=(gx(2)-y1)/(y3-y1); phi(1)=(1-ksi)*(1-eta); phi(2)=(1-ksi)*eta; phi(3)=ksi*eta; phi(4)=ksi*(1-eta);

dxphi(1)=-1/(x2-x1)*(1-eta);

dxphi(2)=-1/(x2-x1)*eta; dxphi(3)=1/(x2-x1)*eta; dxphi(4)=1/(x2-x1)*(1-eta); dyphi(1)=-1/(y3-y1)*(1-ksi); dyphi(2)=1/(y3-y1)*(1-ksi); dyphi(3)=1/(y3-y1)*ksi; dyphi(4)=-1/(y3-y1)*ksi; 3.三角形元的主函数 clear xnum=65; ynum=65; pi=3.141592653;

[point, xn,element,en]=mesh2([0,0],[1,1],xnum,ynum); A=zeros(xn,xn); b=zeros(xn,1);

[gx, w]=gauss(point , element, ei);

[phi, dxphi, dyphi]=shape(point, element, ei, gx); for i=1:3

ii=element(ei,i);

b(ii)=b(ii)+w*8*pi*pi*cos(2*pi*gx(1))*sin(2*pi*gx(2))*phi(i);

for j=1:3

jj=element(ei,j);

A(ii,jj)=A(ii,jj)+w*(dxphi(i)*dxphi(j)+dyphi(i)*dyphi(j));

end end end

if point(i,2)==0||point(i,2)==1 b(i)=0;

A(i,:)=zeros(1,xn); A(i, i)=1; end end for i=1:xn

if point(i,1)==0||point(i,1)==1 b(i)=sin(2*pi*point(i,2)); A(i,:)=zeros(1,xn); A(i, i)=1; end end x=A\b;

%x=inv(A)*b

exact=cos(2*pi*point(:,1)).*sin(2*pi*point(:,2)); %x=inv(A)*b error=0; k=1;

for i=1:xnum for j=1:ynum z(i,j)=x(k); zz(i,j)=exact(k);

error=error+(exact(k)-x(k))^2*1/(xnum-1); k=k+1; end end

fprintf('%f',error); deltax=1.0/(xnum-1);

deltay=1.0/(ynum-1);

mesh([0:deltax:1.0],[0:deltay:1],z) figure;

mesh([0:deltax:1.0],[0:deltay:1],zz) 4.四边形元的主函数 clear xnum=65; ynum=65; pi=3.141592653;

[point, xn,element,en]=mesh1([0,0],[1,1],xnum,ynum); pi=3.141592653; A=zeros(xn,xn); b=zeros(xn,1); for ei=1:en

[gx, w]=gauss1(point , element, ei);

[phi, dxphi, dyphi]=shape1(point, element, ei, gx);

for i=1:4

ii=element(ei,i);

b(ii)=b(ii)+w*8*pi*pi*cos(2*pi*gx(1))*sin(2*pi*gx(2))*phi(i);

for j=1:4

jj=element(ei,j);

A(ii,jj)=A(ii,jj)+w*(dxphi(i)*dxphi(j)+dyphi(i)*dyphi(j));

end end end

for i=1:xn

if point(i,2)==0||point(i,2)==1

b(i)=0;

A(i,:)=zeros(1,xn); A(i, i)=1; end end for i=1:xn

if point(i,1)==0||point(i,1)==1 b(i)=sin(2*pi*point(i,2)); A(i,:)=zeros(1,xn); A(i, i)=1; end end x=A\b;

exact=cos(2*pi*point(:,1)).*sin(2*pi*point(:,2)); %x=inv(A)*b

error=0; k=1;

for i=1:xnum for j=1:ynum z(i,j)=x(k); zz(i,j)=exact(k);

error=error+(exact(k)-x(k))^2*1/(xnum-1); k=k+1; end end

fprintf('%f',error); deltax=1.0/(xnum-1); deltay=1.0/(ynum-1);

mesh([0:deltax:1.0],[0:deltay:1],z) figure;

上机实验（三）

一、实验内容

选用合适的数值方法求解方程：

u

(44x24y2)1u,(x,y)(0,1)(0,1)t

u(0,y,t)u(1,y,t)uy(x,0,t)uy(x,1,t)0;u(x,y,0)sin(x2)cos(y2).

求t0.1,0.5,1.0时，点(

[***********]67

,)(,)(,)(,)(,)、[***********]128128

[***********]3

(,)、(,)、(,)、(,)处的数值解。 [***********]256256

二、实验原理

由于计算机只能存储有限个数据和作有限次运算，所以任何一种用计算机解题的方法，都必须把连续问题（微分方程的边值问题、初值问题）离散化，最终化成有限形式的线性代数方程组。

其求解步骤如下：首先对求解区域进行网格剖分，用有限个网格节点代替连续区域；其次利用微商代替微分的思想，将微分算子离散化，构造逼近微分方程定解问题的差分格式，从而把微分方程的定解问题转化为线性代数方程组问题进行求解。

其求解的微分格式如下：

1nunj,kuj,k

nnnn

unun1j1,k2*uj,kuj1,kj,k12*uj,kuj,k1

*()2244**xj*xj4**yk*ykhh



其中unj,k表示第n层横坐标为j，纵坐标为k的点处的函数值，进行稳定性分析发现，需将区域按256*256来剖分。 三、算法流程

四、计算结果及分析

t取不同值题目中要求的各点的值如下表：

t=0.1s时的图像：

t=0.5时的图像：

t=1时的图像：

分析结果可以发现：随着时间t的增加，所得结果的绝对值减小。

五、总结和讨论

本题中所采用的网格剖分达到了稳定性的要求，算法的时间和空间复杂度都较大，以致计算过程比较复杂。我们需要找到一种方法，既能达到稳定性的要求，又能提高计算效率，这是我们今后努力的方向。

附录：程序源代码

function wentisanchafen(t)

tic

h=1/256;

tao=(1/256)^2;

xnum=1/h+1;

ynum=xnum;

tnum=t/tao;

u0=zeros(xnum,xnum);

u=zeros(xnum,xnum);

x=zeros(1,xnum);

y=zeros(1,xnum);

for j=1:xnum

for n=1:tnum

for j=2:xnum-1

for k=2:ynum-1

x(j)=(j-1)*h;

y(k)=(k-1)*h;

u(j,k)=u0(j,k)+tao*(u0(j+1,k)-2*u0(j,k)+u0(j-1,k)+u0(j,k+1)-2*u0(j,k)+u0(j,k-1))/(4*h*h*(1+pi*x(j)^2+pi*y(k)^2));

u(1,1)=0;

u(1,k)=0;

u(xnum,1)=0;

u(xnum,k)=0;

u(j,1)=u(j,2);

u(j,ynum)=u(j,ynum-1);

end

end

for k=1:ynum

x(j)=(j-1)*h;

y(k)=(k-1)*h;

u0(j,k)=sin(pi*x(j)^2)*cos(pi*y(k)^2); end

end

u0=u;

end

a(1)=u(133,125);

a(2)=u(61,69);

a(3)=u(189,197);

a(4)=u(143,75);

a(5)=u(87,135);

a(6)=u(219,179);

a(7)=u(128,130);

a(8)=u(64,92);

a(9)=u(32,34);

a'

mesh([0:h:1],[0:h:1],u)

toc

2009级数学与应用数学和信息与计算科学专业

偏微分方程数值解上机实验

实验题目 利用有限元方法和有限差分方法求解偏微分方程 完成日期 2012年12月17日 学生姓名 张灵刚 所在班级 1102090 任课教师 王晓东

西北工业大学理学院应用数学系

目录

一．实验目的………………………………………..…….（2） 二．实验要求…………………………………….…...….（2） 三．实验题目……………………………………….………（3） 四．实验二……………………………………………….…（4）

1.实验内容…………………………………….…..（4） 2.实验原理……………………………………….…（4）. 3算法流程………………………………….………（5） 4结果分析…………………………………………..（5）

5总结讨论…………………………………………..（6） 6源程序………………………………………..……（6）

五. 实验三

1.实验内容…………………………………………（17） 2.实验原理…………………………………….….（17）. 3算法流程………………………………………….（18） 4结果分析……………………………………….….（18）.

5总结讨论………………………………………….（21） 6源程序……………………………………………（21）

偏微分方程数值解上机实验报告 实验地点：数学系机房

实验时间：第13—15周，周一、四下午5、6节 实验分数：占期末考试成绩的30% 一、实验目的及意义

掌握有限元方法和有限差分方法的程序实现；学会选择合适的有限差分格式求解一维非线性对流占优的非定常对流扩散问题；学会使用三角线性元和四边形线性元的有限元方法求解二维椭圆方程边值问题，并对计算结果进行收敛性分析；尝试采用有限元方法或有限差分方法实现二维初边值抛物型方程的大规模数值求解。通过实验可以提高学生的动手能力，加深学生对算法的理解。 二、实验要求

在下列给出的三个问题中，最少选择两个问题进行编程实现。要求给出格式的推导过程、算法流程、实现程序、选取的网格参数、以表格或图形的方式给出计算结果、对计算结果进行分析、最后对实验进行总结和讨论。

问题2：用三角线性元和四边形线性元的有限元方法求解方程：

u82cos(2x)sin(2y),(x,y)(0,1)(0,1)u(x,0)u(x,1)0;u(0,y)u(1,y)sin(2y).

取h=1/4, 1/8, 1/16, 1/32, 1/64,比较两种方法的计算精度，并给出数值收敛率.

问题3：选用合适的数值方法求解方程：

u

(44x24y2)1u,(x,y)(0,1)(0,1)t

u(0,y,t)u(1,y,t)uy(x,0,t)uy(x,1,t)0;u(x,y,0)sin(x2)cos(y2).

求t0.1,0.5,1.0时，点(

[***********]67

,)(,)(,)(,)(,)、[***********]128128

[***********]3

(,)、(,)、(,)、(,)处的数值解。 [***********]256256

上机实验（二）

一、实验内容

用三角线性元和四边形线性元的有限元方法求解方程：

u82cos(2x)sin(2y),(x,y)(0,1)(0,1)u(x,0)u(x,1)0;u(0,y)u(1,y)sin(2y).

取h=1/4, 1/8, 1/16, 1/32, 1/64,比较两种方法的计算精度，并给出数值收敛率 二、实验原理

由于计算机只能存储有限个数据和做有限次运算，所以必须把连续问题离散化，有限元法通过把求解偏微分问题转化为变分问题，实质上就是Ritz-Galerkin法，利用有穷维空间近似代替无穷维空间，从而转化为求多元二次函数的极值问题。然后选定单元形状，对求解域进行剖分。构造基函数或单元形状函数，形成有限元空间了，再形成有限元方程，并提供求解有限元方程的有效方法。 三、算法流程

四、计算结果及分析

结果分析：通过给出的上述结果可以发现三种方法相比四边形元的实际求解误差最小，这是由于四边形型元中的基函数中含有

xy项，使得四边形元的误差远比三角形元小，且当空间步长不

断加细时，向前差分的误差逐渐和四边形元接近。而对于收敛率的求解，如果舍去误差，三种方法的收敛率都约为3. 五、总结和讨论

有限元计算的有关问题有：把初值问题化为变分形式，对求解域作网络分割，构造基函数（或单元形状函数），形成有限元

方程。通过本次实验，我懂得了用有限元方法求解初值问题的基本数学思想,理解了线性有限元法的基本原理和方法。另外，我也懂得了按Galerkin方法推导有限元方程的优点，它比Ritz法更加方便直接。我也对虚功原理有了初步的认识。因为Galerkin方法基于虚功原理，所以不但可用于保守场问题，也可使用于非保守场即非驻定问题。 附录：程序源代码 1.三角形元的型函数

function [phi, dxphi, dyphi]=shape(point, element, ei, gx) x1=point(element(ei,1),1); x2=point(element(ei,2),1); x3=point(element(ei,3),1); y1=point(element(ei,1),2); y2=point(element(ei,2),2); y3=point(element(ei,3),2);

S=max(y2-y1,x2-x1)*max(y3-y1,x3-x1);

ksi=((x3*y1-y3*x1)+(y3-y1)*gx(1)+(x1-x3)*gx(2))/S;

eta=((x1*y2-y1*x2)+(y1-y2)*gx(1)+(x2-x1)*gx(2))/S; phi(1)=1-ksi-eta; phi(2)=ksi; phi(3)=eta;

dxphi(1)=-(y3-y1)/S-(y1-y2)/S; dxphi(2)=(y3-y1)/S; dxphi(3)=(y1-y2)/S;

dyphi(1)=-(x1-x3)/S-(x2-x1)/S; dyphi(2)=(x1-x3)/S; dyphi(3)=(x2-x1)/S; 2.四边形元的型函数

function [phi, dxphi, dyphi]=shape1(point, element, ei, gx)

x1=point(element(ei,1),1); x2=point(element(ei,2),1);

x3=point(element(ei,3),1); x4=point(element(ei,4),1); y1=point(element(ei,1),2); y2=point(element(ei,2),2); y3=point(element(ei,3),2); y4=point(element(ei,4),2);

S=max(y2-y1,x2-x1)*max(y4-y1,x4-x1); ksi=(gx(1)-x1)/(x2-x1); eta=(gx(2)-y1)/(y3-y1); phi(1)=(1-ksi)*(1-eta); phi(2)=(1-ksi)*eta; phi(3)=ksi*eta; phi(4)=ksi*(1-eta);

dxphi(1)=-1/(x2-x1)*(1-eta);

dxphi(2)=-1/(x2-x1)*eta; dxphi(3)=1/(x2-x1)*eta; dxphi(4)=1/(x2-x1)*(1-eta); dyphi(1)=-1/(y3-y1)*(1-ksi); dyphi(2)=1/(y3-y1)*(1-ksi); dyphi(3)=1/(y3-y1)*ksi; dyphi(4)=-1/(y3-y1)*ksi; 3.三角形元的主函数 clear xnum=65; ynum=65; pi=3.141592653;

[point, xn,element,en]=mesh2([0,0],[1,1],xnum,ynum); A=zeros(xn,xn); b=zeros(xn,1);

[gx, w]=gauss(point , element, ei);

[phi, dxphi, dyphi]=shape(point, element, ei, gx); for i=1:3

ii=element(ei,i);

b(ii)=b(ii)+w*8*pi*pi*cos(2*pi*gx(1))*sin(2*pi*gx(2))*phi(i);

for j=1:3

jj=element(ei,j);

A(ii,jj)=A(ii,jj)+w*(dxphi(i)*dxphi(j)+dyphi(i)*dyphi(j));

end end end

if point(i,2)==0||point(i,2)==1 b(i)=0;

A(i,:)=zeros(1,xn); A(i, i)=1; end end for i=1:xn

if point(i,1)==0||point(i,1)==1 b(i)=sin(2*pi*point(i,2)); A(i,:)=zeros(1,xn); A(i, i)=1; end end x=A\b;

%x=inv(A)*b

exact=cos(2*pi*point(:,1)).*sin(2*pi*point(:,2)); %x=inv(A)*b error=0; k=1;

for i=1:xnum for j=1:ynum z(i,j)=x(k); zz(i,j)=exact(k);

error=error+(exact(k)-x(k))^2*1/(xnum-1); k=k+1; end end

fprintf('%f',error); deltax=1.0/(xnum-1);

deltay=1.0/(ynum-1);

mesh([0:deltax:1.0],[0:deltay:1],z) figure;

mesh([0:deltax:1.0],[0:deltay:1],zz) 4.四边形元的主函数 clear xnum=65; ynum=65; pi=3.141592653;

[point, xn,element,en]=mesh1([0,0],[1,1],xnum,ynum); pi=3.141592653; A=zeros(xn,xn); b=zeros(xn,1); for ei=1:en

[gx, w]=gauss1(point , element, ei);

[phi, dxphi, dyphi]=shape1(point, element, ei, gx);

for i=1:4

ii=element(ei,i);

b(ii)=b(ii)+w*8*pi*pi*cos(2*pi*gx(1))*sin(2*pi*gx(2))*phi(i);

for j=1:4

jj=element(ei,j);

A(ii,jj)=A(ii,jj)+w*(dxphi(i)*dxphi(j)+dyphi(i)*dyphi(j));

end end end

for i=1:xn

if point(i,2)==0||point(i,2)==1

b(i)=0;

A(i,:)=zeros(1,xn); A(i, i)=1; end end for i=1:xn

if point(i,1)==0||point(i,1)==1 b(i)=sin(2*pi*point(i,2)); A(i,:)=zeros(1,xn); A(i, i)=1; end end x=A\b;

exact=cos(2*pi*point(:,1)).*sin(2*pi*point(:,2)); %x=inv(A)*b

error=0; k=1;

for i=1:xnum for j=1:ynum z(i,j)=x(k); zz(i,j)=exact(k);

error=error+(exact(k)-x(k))^2*1/(xnum-1); k=k+1; end end

fprintf('%f',error); deltax=1.0/(xnum-1); deltay=1.0/(ynum-1);

mesh([0:deltax:1.0],[0:deltay:1],z) figure;

上机实验（三）

一、实验内容

选用合适的数值方法求解方程：

u

(44x24y2)1u,(x,y)(0,1)(0,1)t

u(0,y,t)u(1,y,t)uy(x,0,t)uy(x,1,t)0;u(x,y,0)sin(x2)cos(y2).

求t0.1,0.5,1.0时，点(

[***********]67

,)(,)(,)(,)(,)、[***********]128128

[***********]3

(,)、(,)、(,)、(,)处的数值解。 [***********]256256

二、实验原理

由于计算机只能存储有限个数据和作有限次运算，所以任何一种用计算机解题的方法，都必须把连续问题（微分方程的边值问题、初值问题）离散化，最终化成有限形式的线性代数方程组。

其求解步骤如下：首先对求解区域进行网格剖分，用有限个网格节点代替连续区域；其次利用微商代替微分的思想，将微分算子离散化，构造逼近微分方程定解问题的差分格式，从而把微分方程的定解问题转化为线性代数方程组问题进行求解。

其求解的微分格式如下：

1nunj,kuj,k

nnnn

unun1j1,k2*uj,kuj1,kj,k12*uj,kuj,k1

*()2244**xj*xj4**yk*ykhh



其中unj,k表示第n层横坐标为j，纵坐标为k的点处的函数值，进行稳定性分析发现，需将区域按256*256来剖分。 三、算法流程

四、计算结果及分析

t取不同值题目中要求的各点的值如下表：

t=0.1s时的图像：

t=0.5时的图像：

t=1时的图像：

分析结果可以发现：随着时间t的增加，所得结果的绝对值减小。

五、总结和讨论

本题中所采用的网格剖分达到了稳定性的要求，算法的时间和空间复杂度都较大，以致计算过程比较复杂。我们需要找到一种方法，既能达到稳定性的要求，又能提高计算效率，这是我们今后努力的方向。

附录：程序源代码

function wentisanchafen(t)

tic

h=1/256;

tao=(1/256)^2;

xnum=1/h+1;

ynum=xnum;

tnum=t/tao;

u0=zeros(xnum,xnum);

u=zeros(xnum,xnum);

x=zeros(1,xnum);

y=zeros(1,xnum);

for j=1:xnum

for n=1:tnum

for j=2:xnum-1

for k=2:ynum-1

x(j)=(j-1)*h;

y(k)=(k-1)*h;

u(j,k)=u0(j,k)+tao*(u0(j+1,k)-2*u0(j,k)+u0(j-1,k)+u0(j,k+1)-2*u0(j,k)+u0(j,k-1))/(4*h*h*(1+pi*x(j)^2+pi*y(k)^2));

u(1,1)=0;

u(1,k)=0;

u(xnum,1)=0;

u(xnum,k)=0;

u(j,1)=u(j,2);

u(j,ynum)=u(j,ynum-1);

end

end

for k=1:ynum

x(j)=(j-1)*h;

y(k)=(k-1)*h;

u0(j,k)=sin(pi*x(j)^2)*cos(pi*y(k)^2); end

end

u0=u;

end

a(1)=u(133,125);

a(2)=u(61,69);

a(3)=u(189,197);

a(4)=u(143,75);

a(5)=u(87,135);

a(6)=u(219,179);

a(7)=u(128,130);

a(8)=u(64,92);

a(9)=u(32,34);

a'

mesh([0:h:1],[0:h:1],u)

toc