偏微分方程数值解法

if(aij==NULL)
{
printf("aij fail\n");
aij[i][4]=1.0/(24*N*N)+(-0.5*TSTP)+2*TSTP*1.0/(54*N*N)*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]]);// 1 2
aij[i][5]=1.0/(24*N*N)+(-0.5*TSTP)+2*TSTP*1.0/(54*N*N)*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]]);// 1 3
{
int i;
for(i=1;i<LEE+1;i++)
fe[i][1]=1.0/(18*N*N)*(aryn[a[i][1]]+aryn[a[i][2]]+aryn[a[i][3]])+TSTP*1.0/(54*N*N)*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]])*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]])+TSTP*1.0/(6*N*N)*f((b[a[i][1]].x+b[a[i][2]].x+b[a[i][3]].x)/3.0,(b[a[i][1]].y+b[a[i][2]].y+b[a[i][3]].y)/3.0,(n+1)*TSTP);
{
printf("主元太小");
//主元太小
}
//开始回代
RHS[NG]=RHS[NG]/E[NG][NG];
for(i=NG-1;i>0;i--)
{
double s=0.0;
int j;
for(j=i+1;j<NG+1;j++)
{
s=s+E[i][j]*RHS[j];
}
RHS[i]=(RHS[i]-s)/E[i][i];
{
bd[j]=i;
j++;
}
}Biblioteka Baidu
}
void dealwithbd(double **uk,int bd[]) //近似处理本质边界条件
{
int i;
for(i=bd[1];bd[i]!=0;i++)
uk[bd[i]][bd[i]]=uk[bd[i]][bd[i]]*10e20;
}
void GAUSSELIMINATE(double **E,double RHS[NG+1])//利用高斯消元法解线性方程组
{
if(E[i][k]!=0)
{
double lk=E[i][k]/E[k][k];
int j;
for(j=k+1;j<NG+1;j++)
{
E[i][j]=E[i][j]-lk*E[k][j];
}
RHS[i]=RHS[i]-lk*RHS[k];
}
}
}
if(fabs(E[NG][NG])<1.0e-10)
{
int i;
int k;
for(k=1;k<NG;k++)
{
//选主元
double bmax=0.0;
int ik;
for(i=k;i<NG+1;i++)
{
if(bmax<fabs(E[i][k]))
{
bmax=fabs(E[i][k]);
ik=i;
}
}
if(bmax<1.0e-10)
{
printf("主元太小");
aij[i][6]=1.0/(24*N*N)+2*TSTP*1.0/(54*N*N)*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]]);// 2 3
}
}
void UK(double **uk,int **a,double **aij)//合成总刚度矩阵
{
int i;
aij[i][2]=1.0/(12*N*N)+0.5*TSTP+2*TSTP*1.0/(54*N*N)*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]]);// 2 2
aij[i][3]=1.0/(12*N*N)+0.5*TSTP+2*TSTP*1.0/(54*N*N)*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]]);// 3 3
则 ，有 ， = ，代人（**）即可得到一常微分方程组.
第三步进一步对时间进行离散，得到全离散的逼近格式
对 用差分格式.为此把 等分为n个小区间 ，其长度 , .
这样把求 时刻的近似记为 ， 是 的近似.这里对（**）采用向后的欧拉格式，即
(***)
i=0,1,2…,n-1. =
由于向后欧拉格式为隐式格式且含有非线性项，故相邻两时间步之间采用牛顿迭代，即：
//主元太小
}
//交换第ik行和第k行的元素
if(ik!=k)
{
double t;
for(i=k;i<NG+1;i++)
{
t=E[ik][i];
E[ik][i]=E[k][i];
E[k][i]=t;
}
t=RHS[ik];
RHS[ik]=RHS[k];
RHS[k]=t;
}
//消元
for(i=k+1;i<NG+1;i++)
{if(i%2==1)
{
a[i][1]=(i+1)/2+i/(2*N);
a[i][2]=a[i][1]+1;
a[i][3]=a[i][1]+N+1;
}
if(i%2==0)
{
a[i][3]=i/2+1+(i-1)/(2*N);
a[i][1]=a[i][3]+N+1;
a[i][2]=a[i][3]+N;
L2误差 的阶数为 其中 为时间步长， 为空间步长
这里取N=25,则元素总数LEE = 2*N*N =1250,节点总数NG = (N+1)(N+1) = 676,
取时间步长TSTP = 0.01,时间步数TN = 100,即在t=1时，算得结果为：
即当 = 0.01, =0.04时，
L2误差为0.052853阶数为
3计算各个节点的实际坐标，找出边界点.
4构造单元形状函数，计算雅可比矩阵并利用高斯求积公式计算单元刚度矩阵和单元列阵.
5合成总刚度矩阵和总列阵.
6处理本质边界条件，求解线性方程组.
7在第n个时间步利用牛顿迭代法算得第n+1个时间步的数值解，取 .
8计算各个时间步的有限元数值解和L2误差.
四、误差分析
b[i].y=(i/(N+1))*(1.0/N);
}
}
void boundnote(int bd[],struct xy b[])//找出边界点
{
int i,j=1;
for(i=1;i<NG+1;i++)
{
if(b[i].x==0||b[i].x==1.0||b[i].y==0||b[i].y==1.0)
uk[a[i][1]][a[i][3]]+=aij[i][5];
uk[a[i][3]][a[i][1]]+=aij[i][5];
uk[a[i][2]][a[i][3]]+=aij[i][6];
uk[a[i][3]][a[i][2]]+=aij[i][6];
}
}
void FE(double **fe,double aryn[],double aryk[],int n,int **a,struct xy b[]) //计算单元列阵
return z;
}
double u(double x,double y,double t)//方程的准确解
{
double z;
z=100*(t+1)*x*y*(x-1)*(y-1);
return z;
}
void II(int **a) //节点的局部编码与总体编码
{
int i;
for(i=1;i<LEE+1;i++)
一、问题
用有限元方法求下面方程的数值解
in
on
in
二、问题分析
第一步利用Green公式，求出方程的变分形式
变分形式为：求 ，使得
（*）
以及 .
第二步对空间进行离散，得出半离散格式
对区域 进行剖分，构造节点基函数，得出有限元子空间： ，则（*）的Galerkin逼近为：
，求 ，使得
（**）
以及 ， 为初始条件 在 中的逼近，设 为 在 中的插值.
{
int i;
double *aryk;
aryk=(double *)calloc(NG+1,sizeof(double));
for(i=1;i<NG+1;i++)
aryk[i]=aryn[i];
double boot;
double **uk;
double *ufe;
double **aij;
double **fe;
#define TSTP 0.01 //时间步长
#define TN 100 //时间迭代步数
#define J 1.0/(N*N) //雅可比行列式的绝对值
double u0(double x,double y) //初值函数u0
{
double z;
z=100*x*y*(x-1)*(y-1);
return z;
if(uk[i]==NULL)
{
printf("uki fail\n");
exit(1);
}
}
ufe=(double *)calloc(NG+1,sizeof(double));
if(ufe==NULL)
{
printf("ufe fail\n");
exit(1);
}
aij=(double **)calloc(LEE+1,sizeof(double *));
其中用 作为牛顿迭代的初值.
三、计算流程
取 作为方程的准确解，算得初值函数：
右端函数：
+
.
取 .
1对 作三角剖分：分别沿x , y方向对[0,1]区间作N等分并将每个小正方形沿对角线分为两个三角形.
2对节点进行总体编码（这里按先沿x方向再沿y方向进行顺序编码），并建立局部编码和总体编码之间的对应关系.
}
}
}
struct xy {double x;double y;};
void XY(struct xy b[]) //节点实际坐标
{
int i;
for(i=1;i<NG+1;i++)
if(i%(N+1)==0)
{
b[i].x=1;
b[i].y=(i/(N+1)-1)*(1.0/N);
}
else
{
b[i].x=(i%(N+1)-1)*(1.0/N);
do
{
uk=(double **)calloc(NG+1,sizeof(double *));
if(uk==NULL) //判断内存分配是否成功
{
printf("uk fail\n");
exit(1);
}
for(i=0;i<NG+1;i++)
{
uk[i]=(double *)calloc(NG+1,sizeof(double));
附C程序
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#define N 25 //沿x y方向的等分数
#define ND 3 //一个元素的节点个数
#define LEE 2*N*N //元素总数
#define NG (N+1)*(N+1) //节点总数
}
void UFE(double ufe[],double **fe,int **a)//合成总列阵
{
int i;
for(i=1;i<LEE+1;i++)
{
ufe[a[i][1]]+=fe[i][1];
ufe[a[i][2]]+=fe[i][1];
ufe[a[i][3]]+=fe[i][1];
}
}
void NEWTONITERATIVE(double aryn1[],double aryn[],int n,int **a,struct xy b[],int bd[])//牛顿迭代
}
}
void AIJ(double **aij,double aryk[],int **a) //计算单元刚度矩阵
{
int i;
for(i=1;i<LEE+1;i++)
{
aij[i][1]=1.0/(12*N*N)+TSTP+2*TSTP*1.0/(54*N*N)*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]]);// 1 1
}
double f(double x,double y,double t)//右端函数f
{
double z;
z=100*x*y*(x-1)*(y-1)-200*(t+1)*y*(y-1)-200*(t+1)*x*(x-1)+10000*(t+1)*(t+1)*x*x*y*y*(x-1)*(x-1)*(y-1)*(y-1);
for(i=1;i<LEE+1;i++)
{
uk[a[i][1]][a[i][1]]+=aij[i][1];
uk[a[i][2]][a[i][2]]+=aij[i][2];
uk[a[i][3]][a[i][3]]+=aij[i][3];
uk[a[i][1]][a[i][2]]+=aij[i][4];
uk[a[i][2]][a[i][1]]+=aij[i][4];