有限元方法求解边值问题

1、有限元方法求解边值问题一、问题用有限元方法求解边值问题-u + u = (1 + 7t2)sin(开 x),0 < % < 1,u(0) = 0, u(l) = 0 .二.求解过程已知该问题的精确解为u(x) = sin( n x).有限元方法：1、单元剖分 将区间0,1作n+1等分，记h = 1/5 + 1).2、构造有限元空间构造卅(0,1)的有限维子空间.年严，< x <3、有限元空间的基函数0心)=：宀,xi<x< xi+1 i = 1, 2,0,其他4、有限元方程组/q(01，01)q(0i，02)q(02, 01)q(02, 02)q(02, 0

2、3) a(0n-l，0n_2)a(0n-l> 0n-l) q(0n，0归)/ (a 01) (a 02)(a 0-1) (f，0n) / 由此方程组解出sc?,心为问题的有限元解.q(0n-八 n-1q(00j / cn /三.结果当n二1时，基函数为(2. 0 < % < | ,01 (x) = 2(2(1 %) , < x <1.有限元方程组为解得有限元解为& 6667c =8.8106 ,5 =1.0166 .ui(x) =1. 01660(x).当n二4时，基函数为(5%,0 < x < | ,0i(x) = « 5(|-兀)

3、 |<x<| 0 , x不属于0,- 'l 5l-x<l，02(x)= « 瓷-%), |<x<| '、0 , x不属于£,| ,5（兀-|），lx<l，0,x） = 5qx）' ”杲2 4'l5,5（rf 3、3,45（光一»5-x<i404（x） =5（1 -x）z - <x <1 f0 ,尢不属于i，1 .5有限元方程组为/ 50.6667i -24.8333-24.833350.6667-24.8333-24.833350.6667-24.8333-24.833350.6

4、6676.1816e + 001.0002e + 011.0002e + 016.1816e + 00解得c = 0.58953, c2 = 0.95389, c3 = 0.95389, c4 = 0.58953. 有限元解为u4w = sf=iq0i（x）.当n二8时，基函数为9%,0 < % < | ,0i（x）=9（|-x）, i<x<| ,、0 ,兀不属于0,| ,%-|）,|<x<| ,03（x） = * 9（才一%）鳥 sx 今,0 ,尢不属于篙，'9（x -扌）注 * < £ 05（x） = < 9（彳一%）鳥 s

5、xsf, 0 ,尢不属于冷月,、l9 99（% -） ,- < x < -, 9八 9907（x） = « 9（善一x） , 1<x<1 ,0,兀不属于e月,l9 9|<x<|,02（x）= 9（|-x）, |<x<| , 、0 '兀不属于耳，i'04（x）= < 9（|_x）-x-|， 0 ,%不属于|，|，”9（9， |3v 齐06（x） = « g-x） , |<<5， 0 , %不属于|冷，%-9，討龙v齐o08（x） =9（1 -x）, - <x< 1 ,0 ,兀不属于|

6、,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.9162 19.318210.59610.5969.3182 (6.9162、36800/解得c = 0.34234, c2 = 0.64339, c3 = 0.86683, c4 = 0.98572, c5 = 0.98572, c6 = 0.8

7、6683, c7 = 0.64339, c8 = 0.34234. 有限元解为8l=s=8c；0i(x)團6 6有限元解曲线u4(x)和u8(x)四、程序求系数clc, clear %有限元方程组 ax 二 bn 二 input. ('plaose input n:'); h=l/(n+1);x=0+(0:n+l)*h;ux=sin(pi*x(2:n+l) % 精确 解a=a/hfor i二2:n+1fl=sym(，(l+pi"2)*sin(pi*x)');f2=sym(, (1+p厂2)*sin(pi*x)*x')fl=sym(，l/h+h*x 2+

8、1/h+h*(lx)al=eval (int (f 1, 0, 1);f2二sym('-l/h+h*(l-x)*x');a2=eval (int (f2, 0, 1);a=；for i=l:na(i, i)=a.l;endfor i=l:n-la(i+l, i)=a2;a(i, i+l)=a2;endpl=eval (int (f2, x (it), x (it) +h);p2=eval (int (fl, x (i-1), x(il)+h);p3=eval (int (fl, x(i), x(i)+h);p4=eval (int (f2, x (i), x (i) +h)；b

9、(i-1)=l/h*pl-x(i-1)/h*p2+(l+x(i)/h)*p3-l/h*p4;endb二b/h;b 二 b'c=inv(a) *b %求有限元系数基函数文件function f=fpp(x, i, n) % righth=l/(n+l);y=0+(0:n+l)*h;if x>=y(i) & x<y(i+l)f 二(x-y(i)/h;else if x>=y(i+1) & x<=y(i+2)f 二(y (i+2) - x)/h;else f=0;endend图6. 7clc, clear %有限元方程组ax二bfor k二1:2n 二

10、 input ('please in put n:'); h=l/(n+l);x=0+(0:n+l)*h;fl二sym(' l/h+h*x"2+l/h+h*(l-x厂2'); al二eval (int (fl, 0, 1);f2=sym(，-l/h+h*(l-x)*x');a2=eval (int (f2, 0, 1);a=;for i=l:na(i, i)=al;endfor i二l:nta(i + l, i)=a2;a(i, i+l)=a2;enda=a/h;b=;for i=2:n+lfl=sym(，(1+p厂2)*sin(pi*x)&#

11、39;);f2=sym(，(l+pi"2)*sin(pi*x)*x);p2=eval(int(fl, x(il), x(il)+h); p3=eval (int (f 1, 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+ (1+x (i) /h )*p3t/h*p4;endb=b/h;b二b'c二inv(a)*b;%求有限元系数ux=;x=0:0. 01:1;ux二sin(pi*x);% 精确解%求有限元解p=；for j=l:nfor i=l:length(x)p(i, j)=fpp(x(i), j,n);endendp二p' ;u=c *p;e=abs (u-ux);if n=4plot (x, e,'；hold onen