做网站策划计划书,正规东莞网站建设,上海网站建设服务宁德,网站建设是自己做好还是外包目录
一、研究目标
二、理论推导
三、算例实现
四、结论 一、研究目标 我们在前几节中介绍了Poisson方程的边值问题#xff0c;接下来对椭圆型偏微分方程的混合边值问题进行探讨#xff0c;研究对象为#xff1a; 其中#xff0c;为矩形区域#xff0c;为上的连续函数…目录
一、研究目标
二、理论推导
三、算例实现
四、结论 一、研究目标 我们在前几节中介绍了Poisson方程的边值问题接下来对椭圆型偏微分方程的混合边值问题进行探讨研究对象为 其中为矩形区域为上的连续函数是上的单位法向量从而表示方向导数为的连续函数且非负。对于矩形区域而言其边界上的法向量没有统一的表达式需要对四条边界线段分别讨论。可见 在左边界边界条件为 在右边界边界条件为 在下边界边界条件为 在上边界边界条件为。 于是公式1可以具体写为 二、理论推导 首先进行矩形区域的等距剖分得到各网格节点且。然后弱化原方程使其仅在离散节点上成立从而有 将上式中的一阶、二阶偏导数分别用关于一阶导数的中心差商和关于二阶导数的中心差商来近似得 用数值解代替精确解并忽略高阶项得到以下数值格式 其中。该格式的局部截断误差为同时需要处理其中的下标越界问题。由于函数的连续性我们认为公式2中的第1式对i0也成立即有 再成公式2的第2式中解出代入上式得 同样设公式2中的第1式分别对成立再分别从公式2的第3、4、5式中解出代入前面刚得到的方程就可以处理掉越界下标得到以下格式 至此我们共有m1xn1个待求量,而现有m-1xn-1个关于内节点的方程2n-1个关于左、右边界上的节点不含端点的方程及2m-1个关于上、下边界上的节点也不含端点的方程还需要补充 个方程也就是关于矩形区域的4个顶点的方程。为此设公式2中第1式对成立即 然后再从公式2的第2式和第4式中分别解出和代入公式3中得到左下顶点处的方程 上式可以简化为 同样设公式2中第1式分别对、和成立然后再从公式2中第3式和第4式解出分别代入刚才得到的3个方程就得到的右下顶点、左上顶点、和右上顶点处的方程这样我们就有了完整的处理带导数边界条件的椭圆型方程的数值格式 记有 为计算分别可将上述方程组写成矩阵形式。 首先将公式4中第4、6、7式写成 上面的式子可以简记为 其中为m1阶单位矩阵且C为公式5最左端的三对角矩阵为公式5右端的向量。接着公式4中的第1,2,3式可以写成
, 该式可以简记为 其中为公式6中的三对角矩阵为公式6右端的向量。最后公式4的第5、8、9式可以写成 该式可以简记为 其中D为公式7中的三对角矩阵为公式7右端的向量。于是由公式5、6、7可得到公式4写成块三对角矩阵为 采用Gauss-Seidel迭代法求解公式8。
三、算例实现 用差分格式-公式8求解椭圆型方程混合边值问题 已知精确解为。分别取步长为和输出6个节点和处的数值解并给出误差要求在各节点处最大误差的迭代误差限为。
代码如下 #include cmath
#include stdlib.h
#include stdio.h
#define pi 3.14159265359int main(int argc, char*argv[])
{int m,n,i,j,k;double xa,xb,ya,yb,dx,dy,alpha,beta,gamma,maxerr;double *x,*y,**u,**v,**lambda,*d,kexi,eta,temp;double f(double x, double y);double lambda_function(double x, double y);double phi1(double y);double phi2(double y);double psi1(double x);double psi2(double x);double exact(double x, double y);xa0.0;xb2.0;ya0.0;yb1.0;m64;n32;printf(m%d,n%d.\n,m,n);dx(xb-xa)/m;dy(yb-ya)/n;beta1.0/(dx*dx);gamma1.0/(dy*dy);alpha2*(betagamma);kexi2.0/dx;eta2.0/dy;x(double*)malloc(sizeof(double)*(m1));for(i0;im;i)x[i]xai*dx;y(double*)malloc(sizeof(double)*(n1));for(j0;jn;j)y[j]yaj*dy;u(double**)malloc(sizeof(double*)*(m1));v(double**)malloc(sizeof(double*)*(m1));lambda(double**)malloc(sizeof(double*)*(m1));for(i0;im;i){u[i](double*)malloc(sizeof(double)*(n1));v[i](double*)malloc(sizeof(double)*(n1));lambda[i](double*)malloc(sizeof(double)*(n1));}for(i0;im;i){for(j0;jn;j){u[i][j]0.0;v[i][j]0.0;lambda[i][j]lambda_function(x[i],y[j]);}}d(double*)malloc(sizeof(double)*(m1));k0;do{maxerr0.0;for(i0;im;i)d[i]f(x[i],y[0])-eta*psi1(x[i]);d[0]d[0]-kexi*phi1(y[0]);d[m]d[m]kexi*phi2(y[0]);v[0][0](d[0]2*gamma*u[0][1]2*beta*u[1][0])/(alpha(kexieta)*lambda[0][0]);for(i1;im;i)v[i][0](d[i]2*gamma*u[i][1]beta*(v[i-1][0]u[i1][0]))/(alphaeta*lambda[i][0]);v[m][0](d[m]2*gamma*u[m][1]2*beta*v[m-1][0])/(alpha(kexieta)*lambda[m][0]);for(j1;jn;j){for(i0;im;i){d[i]f(x[i],y[j]);}d[0]d[0]-kexi*phi1(y[j]);d[m]d[m]kexi*phi2(y[j]);v[0][j](d[0]gamma*(u[0][j1]v[0][j-1])2*beta*u[1][j])/(alphakexi*lambda[0][j]);for(i1;im;i){v[i][j](d[i]gamma*(v[i][j-1]u[i][j1])beta*(v[i-1][j]u[i1][j]))/alpha;}v[m][j](d[m]gamma*(v[m][j-1]u[m][j1])2*beta*v[m-1][j])/(alphakexi*lambda[m][j]);}for(i0;im;i)d[i]f(x[i],y[n])eta*psi2(x[i]);d[0]d[0]-kexi*phi1(y[n]);d[m]d[m]kexi*phi2(y[n]);v[0][n](d[0]2*gamma*v[0][n-1]2*beta*u[1][n])/(alpha(kexieta)*lambda[0][n]);for(i1;im;i)v[i][n](d[i]2*gamma*v[i][n-1]beta*(v[i-1][n]u[i1][n]))/(alphaeta*lambda[i][n]);v[m][n](d[m]2*gamma*v[m][n-1]2*beta*v[m-1][n])/(alpha(kexieta)*lambda[m][n]);for(i0;im;i){for(j0;jn;j){tempfabs(u[i][j]-v[i][j]);if(tempmaxerr)maxerrtemp;u[i][j]v[i][j];}}kk1;}while((maxerr0.5*1e-10)(k1e8));printf(k%d\n,k);km/4;for(ik;im;iik){printf((%.2f,0.25), y%f, err%.4e.\n,x[i],u[i][n/4],fabs(exact(x[i],y[n/4])-u[i][n/4]));}km/4;for(ik;im;iik){printf((%.2f,0.50), y%f, err%.4e.\n,x[i],u[i][n/2],fabs(exact(x[i],y[n/2])-u[i][n/2]));}for(i0;im;i){free(u[i]);free(v[i]);free(lambda[i]);}free(u);free(v);free(lambda);free(x);free(y);free(d);return 0;
}double f(double x, double y)
{return (pi*pi-1)*exp(x)*sin(pi*y);
}
double lambda_function(double x, double y)
{return x*xy*y;
}
double phi1(double y)
{return (1-y*y)*sin(pi*y);
}
double phi2(double y)
{return (5.0y*y)*exp(2.0)*sin(pi*y);
}
double psi1(double x)
{return pi*exp(x);
}
double psi2(double x)
{return -pi*exp(x);
}
double exact(double x, double y)
{return exp(x)*sin(pi*y);
} 当时计算结果如下
m32,n16.
k4323
(0.50,0.25), y1.152179, err1.3643e-02.
(1.00,0.25), y1.911016, err1.1100e-02.
(1.50,0.25), y3.162159, err6.8738e-03.
(0.50,0.50), y1.638607, err1.0115e-02.
(1.00,0.50), y2.711255, err7.0265e-03.
(1.50,0.50), y4.479936, err1.7526e-03.
当时计算结果如下
m64,n32.
k16007
(0.50,0.25), y1.162412, err3.4097e-03.
(1.00,0.25), y1.919341, err2.7745e-03.
(1.50,0.25), y3.167313, err1.7193e-03.
(0.50,0.50), y1.646193, err2.5286e-03.
(1.00,0.50), y2.716524, err1.7575e-03.
(1.50,0.50), y4.481249, err4.3972e-04.
四、结论 从计算结果可知当步长减小为1/2时误差减小为原来的1/4可见该数值格式是二阶收敛的。
