偏微分方程數(shù)值解法

一、問(wèn)題ffinoninuuu2tu0uxu0二、問(wèn)題分析第一步利用Green公式，求出方程的變分形式變分形式為：求uL20,T;H1，使得0u,vu,vu2,vf,v vH1 （）t 0和 uxu.0第二步對(duì)空間進(jìn)行離散，得出半離散格式對(duì)區(qū)域

,

,,

，那么（*）的Galerkin逼近為：

h 1 2 NGt0,T，求ut,xVh

H1，使得0d,vdt h h

t,vut2,vf,vh h h

Vh h

（**）和u0u

，u0,

為初始條件u0

在V中的逼近，設(shè)uh

為u0,h

在V.h則t0，有u

tNGt，u

=NG

，代人（**）即可取得一常微分方程組.h i1

0,

i1

0i i第三步進(jìn)一步對(duì)時(shí)刻進(jìn)行離散，取得全離散的逼近格式對(duì)du等分為n個(gè)小區(qū)間,tdt i1

，其長(zhǎng)度tti

ti1

T,tn

T. 如此把求t時(shí)刻的近似記為ui

的近似.那個(gè)地址對(duì)（**）采納向后的歐拉格式，即i1

h 0 ui1ui,v

ui1,v ui1 ,

f

,v V

(***)t h h h

h h h h

ih h hi=0,1,2…,n-1. u0=h

0,h由于向后歐拉格式為隱式格式且含有非線性項(xiàng)，故相鄰兩時(shí)刻步之間采納牛頓迭代，即：u un1 k h ,

,v

2u

u,

2,

f,v h k h

kk h

kh n hu其頂用un1作為牛頓迭代的初值.hu三、計(jì)算流程取u100t1xyx1y1作為方程的準(zhǔn)確解，算得初值函數(shù)：u 100xyx1y10右端函數(shù)：f100xyx1y1200t1xx1200t1yy1+10000t2x2y2x2y2.取y;0x1,0y對(duì)x,y方向?qū)0,1]N等分并將每一個(gè)小正方形沿對(duì)角線分為兩個(gè)三角形.對(duì)節(jié)點(diǎn)進(jìn)行整體編碼（那個(gè)地址按先沿x方向再沿y方向進(jìn)行順序編碼..陣...n個(gè)時(shí)刻步利用牛頓迭代法算得第n+1個(gè)時(shí)刻步的數(shù)值解，取u0h

u .0,h計(jì)算各個(gè)時(shí)刻步的有限元數(shù)值解和L2誤差. L2誤差utn

unh0

的階數(shù)為th2

其中t為時(shí)刻步長(zhǎng)，h為空間步長(zhǎng)那個(gè)地址取N=25,那么元素總數(shù)LEE=2*N*N=1250,節(jié)點(diǎn)總數(shù)NG=(N+1)(N+1)=676,取時(shí)刻步長(zhǎng)TSTP=,時(shí)刻步數(shù)TN=100,即在t=1時(shí)，算得結(jié)果為：即當(dāng)t=,h=時(shí)，L2誤差為階數(shù)為th2附C程序#include""#include""#include""#defineN25 =1;b[i].y=(i/(N+1)-1)*N);}else{}}

b[i].x=(i%(N+1)-1)*N);b[i].y=(i/(N+1))*N);voidboundnote(intbd[],structxyb[])==0||b[i].x==||b[i].y==0||b[i].y=={bd[j]=i;j++;}}}void dealwithbd(double **uk,int bd[])+b[a[i][2]].x+b[a[i][3]].x)/,(b[a[i][1]].y+b[a[i][2]].y+b[a[i][3]].y)/,(n+1)*TSTP);}voidUFE(doubleufe[],double**fe,int**a),b[i].y);double*aryn1;for(n=0;n<TN;n++){aryn1=(double*)calloc(NG+1,sizeof(double));NEWTONITERATIVE(aryn1,aryn,n,a,b,bd);for(i=0;i<NG+1;i++)aryn[i]=aryn1[i];free(aryn1);}for(i=1;i<NG+1;i++)fesolution[i]=aryn[i];printf("\n第%d步的結(jié)果為：\n\n",n);printf("節(jié)點(diǎn)整體編碼 有限元值 精準(zhǔn)值 誤差 for(j=N/2;j<NG+1;j+=2*(N+1))printf("%8d%15f%15f%15f\n",j,fesolution[j],u(b[j].x,b[j].y,TSTP*TN),fesolution[j]-u(b[j].x,b[j].y,TSTP*TN