有限元方法求解边值问题

《有限元方法求解边值问题》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

1、有限元方法求解边值问题一、问题用有限元方法求解边值问题-u+u=(1+7T2)sin(开x),0<%<1,u(0)=0,u(l)=0.二.求解过程已知该问题的精确解为u(x)=sin(nx).有限元方法：1、单元剖分将区间［0,1］作n+1等分，记h=1/5+1).2、构造有限元空间构造卅(0,1)的有限维子空间.年严，

2、_2)a(0n-l>0n-l)Q(0n，0归)/(A01)(A02)■•(A0-1)(f，0n)/由此方程组解出SC?,…心为问题的有限元解.Q(0n-八^n-1Q(0〃0J/Cn/三.结果当n二1时，基函数为(2.0<%<

3、,01(x)=]]2(2(1—%),—

4、,0i(x)=«5(

5、-兀)‘

6、

7、‘0,x不属于0,-‘'L5」l-x

8、

9、

10、'、0,x不属于£,

11、,5（兀-

12、），l~x

13、.58953.有限元解为u4W=Sf=iQ0i（x）.当n二8时，基函数为9%,0<%<

14、,0i（x）=9（

15、-x）,i

16、,、0,兀不属于0,

17、,%-

18、）,

19、

20、,03（x）=*9（才一%）鳥Sx今,0,尢不属于篙，'9（x-扌）’注*<£‘05（x）=<9（彳一%）鳥Sxsf,0,尢不属于［冷月,、L99」9（%—-）,-

21、

22、,02（x）=9（

23、-X）,

24、

25、,、0'兀不属于耳，I]'04（X）=<9（

26、_X）^-X-

27、，0,%

28、不属于[

29、，

30、]，”9（—9，

31、3V齐06（x）=«^g-x）,

32、<^<5，0,%不属于[

33、冷，%-9，討龙V齐o08（x）=9（1-x）,-

34、,1]・有限元方程组为/162.667-80.833-80.833162.667-80.833-80.833162.667-80.833-80.833162.667-80.833-80.833162.667-80.833-80.833162.667-80.833-80.833162.667-80.833-80.833162.667C5/3.68006.916219.31

35、8210.59610.5969.3182(6.9162、3・6800/解得C]=0.34234,c2=0.64339,c3=0.86683,c4=0.98572,c5=0.98572,c6=0.86683,c7=0.64339,c8=0.34234.有限元解为8l=s=8c；0i(x)・團66有限元解曲线u4(x)和u8(x)四、程序求系数clc,clear%有限元方程组Ax二bn二input.('Plaoseinputn:');h=l/(n+1);x=0+(0:n+l)*h;ux=sin(pi*x(2:n+l))%精确解a=a/hf

36、ori二2:n+1fl=sym(，(l+pi"2)*sin(pi*x)');f2=sym(,((1+p厂2)*sin(pi*x))*x')fl=sym(，l/h+h*x2+1/h+h*(l~x)al=eval(int(f1,0,1));f2二sym('-l/h+h*(l-x)*x');a2=eval(int(f2,0,1));a=[]；fori=l:na(i,i)=a.l;endfori=l:n-la(i+l,i)=a2;a(i,i+l)=a2;endpl=eval(int(f2,x(iT),x(iT)+h));p2=eval(int

37、(fl,x(i-1),x(i~l)+h));p3=eval(int(fl,x(i),x(i)+h));p4=eval(int(f2,x(i),x(i)+h))；b(i-1)=l/h*pl-x(i-1)/h*p2+(l+x