资源描述:
《有限元方法求解边值问题》由会员上传分享,免费在线阅读,更多相关内容在工程资料-天天文库。
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-x8、
9、10、'、0,x不属于£,
11、,5(兀-
12、),l~x13、.58953.有限元解为u4W=Sf=iQ0i(x).当n二8时,基函数为9%,0<%<
14、,0i(x)=9(
15、-x),i16、,、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