有限元方法求解初边值问题

1 9
1
9(������ − 9) , ∅2 (x) = 9 (9 − ������) , { 0 ,
3 5 3
1
1 9 2 9
≤ ������ < 9 , ≤ ������ ≤ 9 ,
1 3 9 9 4 3
2
∅1 (x) = 9 (9 − ������) , { 0 ,
4 2
≤ ������ ≤ 9 ,
≤ ������ <
3 5 4
, ∅4 (x) =
5(������ − ) ,
5
3
3 5 4 5
≤ ������ <
3
4 5
,
∅3 (x) =
5 (5 − ������) , { 0 ,
4
≤ ������ ≤ 5 ,
2 4
5(1 − ������) , { 0 ,
≤ ������ ≤ 1 ,
������不属于 [5 , 5] ,
6 8
≤ ������ ≤ 1 ,
������不属于 [9 , 9] ,
������不属于 [ , 1] .
有限元方程组为 162.667 −80.833 −80.833 162.667 −80.833
(
������1 ������ −80.833 2 ������3 162.667 −80.833 ������4 −80.833 162.667 −80.833 ������5 −80.833 162.667 −80.833 −80.833 ������6 −80.833 162.667 ������7 −80.833 162.667 −80.833 −80.833 162.667 ) (������8 )
9(������ − 9) , ∅4 (x) = 9 (9 − ������) , { 0 ,
5 9
≤ ������ < 9 , ≤ ������ ≤ 9 ,
3 5 6 9 7 9 5
������不属于 [9 , 9] ,
4 9 5 9
������不属于 [9 , 9] ,
5 9 6 9
9(������ − ) , ∅5 (x) = 9 ( − ������) , { 0 ,
a=a/h for i=2:n+1 f1=sym('(1+pi^2)*sin(pi*x)'); f2=sym('((1+pi^2)*sin(pi*x))*x') ; p1=eval(int(f2,x(i-1),x(i-1)+h)); p2=eval(int(f1,x(i-1),x(i-1)+h)); p3=eval(int(f1,x(i),x(i)+h)); p4=eval(int(f2,x(i),x(i)+ h)); b(i-1)=1/h*p1-x(i-1)/h*p2 +(1+x(i)/h)*p3-1/h*p4; end b=b/h;
ℎ������ ������������+1 −������ ℎ������+1
, ������������−1 ≤ ������ < ������������ , ������������ ≤ ������ ≤ ������������+1 0, 其他 ������1 ������2 ⋮ i = 1，2， ⋯ ，n.
2
������1 6.1816e + 00 ������2 −24.8333 ) (������ )=(1.0002e + 01). 3 1.0002e + 01 50.6667 −24.8333 ������4 6.1816e + 00 −24.8333 50.6667
0 ≤ ������ < 9 ,
四、程序 求系数 clc,clear % 有限元方程组 Ax=b n=input('Plaese input n: '); h=1/(n+1); x=0+(0:n+1)*h; ux=sin(pi*x(2:n+1)) % 精确 解 f1=sym('1/h+h*x^2+1/h+h*(1-x) ^2'; a1=eval(int(f1,0,1)); f2=sym('-1/h+h*(1-x)*x'); a2=eval(int(f2,0,1)); a=[]; for i=1:n a(i,i)=a1; end for i=1:n-1 a(i+1,i)=a2; a(i,i+1)=a2; end
1 ≤ ������ ≤ 1 . 2
1 , 2
8.6667������1 =8.8106 , 解得 ������1 =1.0166 . 有限元解为 ������1 (������) =1.0166∅1 (x). 当 n=4 时，基函数为 5������ ,
2
0 ≤ ������ < 5 ,
1 5
1
5(������ − 5) , ∅2 (x) = 5 (5 − ������) , { 0 ,
3
1
1 5 2 5
≤ ������ < 5 , ≤ ������ ≤ 5 , Nhomakorabea1 3 3
2
∅1 (x) = 5 (5 − ������) , { 0 ,
≤ ������ ≤ 5 ,
2
2
������不属于 [0, 5] ,
������不属于 [5 , 5] ,
5(������ − ) ,
5
2
2 5 3 5
2 9 2 3
2
������不属于 [0, ] ,
9 3 9
������不属于 [ , ] ,
3 9 4 9
9(������ − 9) , ∅3 (x) = 9 (9 − ������) , { 0 ,
4 9 6 9
≤ ������ < 9 , ≤ ������ ≤ 9 ,
2 4 5 9 6 9 4
⋱ ������������−1 ������(∅������−1 , ∅������−2 ) ������(∅������−1 , ∅������−1 ) ������(∅������−1 , ∅������ ) ������(∅������ , ∅������−1 ) ������(∅������ , ∅������ ) ) ( ������������ )
有限元方法求解初边值问题
一、问题 用有限元方法求解边值问题 −������′′ + ������ = (1 + ������ 2 )sin(πx), 0 < ������ < 1 , ������(0) = 0, ������(1) = 0 . 二． 求解过程 已知该问题的精确解为 ������(������) = sin(πx) . 有限元方法： 1、单元剖分 将区间 [0,1] 作 n+1 等分，记 h = 1/(������ + 1). 1 2、构造有限元空间 构造������0 (0,1)的有限维子空间.
������−������������−1
3、有限元空间的基函数 ∅������ (������) = { 4、有限元方程组 ������(∅1 , ∅1 ) ������(∅1 , ∅2 ) ������(∅2 , ∅1 ) ������(∅2 , ∅2 ) ������(∅2 , ∅3 ) ⋱ (
b=b' c=inv(a)*b % 求有限元系数 p2=eval(int(f1,x(i-1),x(i-1)+h)); p3=eval(int(f1,x(i),x(i)+h)); p4=eval(int(f2,x(i),x(i)+h)); b(i-1)=1/h*p1-x(i-1)/h*p2+(1+x(i)/h )*p3-1/h*p4; end b=b/h; b=b'; c=inv(a)*b; % 求有限元系数 ux=[]; x=0:0.01:1; ux=sin(pi*x); % 精确解 % 求有限元解 P=[]; for j=1:n for i=1:length(x) P(i,j)=fpp(x(i),j,n); end end P=P';u=c'*P; E=abs(u-ux); if n==4 plot(x,E,'-'); hold on end if n==8 plot(x,E,'-x'); end legend('|u(x)-u4(x)|','|u(x)-u8(x)| '); title('图 6.7 有限元解误差 |u(x)-un(x)|曲线'); k=k+1; end
������不属于 [5 , 1] .
有限元方程组为 50.6667 −24.8333 −24.8333 50.6667 ( −24.8333 解得 ������1 = 0.58953, ������2 = 0.95389, ������3 = 0.95389, ������4 = 0.58953. 有限元解为 ������4 (������) = ∑4 ������=1 ������������ ∅������ (x). 当 n=8 时，基函数为 9������ ,
基函数文件 function f=fpp(x,i,n) % right h=1/(n+1); y=0+(0:n+1)*h; if x>=y(i) & x<y(i+1) f=(x-y(i))/h; else if x>=y(i+1) & x<=y(i+2) f=(y(i+2)-x)/h; else f=0; end end 图 6.7 clc,clear % 有限元方程组 Ax=b for k=1:2 n=input('Please input n: '); h=1/(n+1); x=0+(0:n+1)*h; f1=sym('1/h+h*x^2+1/h+h*(1-x)^2'); a1=eval(int(f1,0,1)); f2=sym('-1/h+h*(1-x)*x'); a2=eval(int(f2,0,1)); a=[]; for i=1:n a(i,i)=a1; end for i=1:n-1 a(i+1,i)=a2; a(i,i+1)=a2; end a=a/h; b=[]; for i=2:n+1 f1=sym('(1+pi^2)*sin(pi*x)'); f2=sym('((1+pi^2)*sin(pi*x))*x'); p1=eval(int(f2,x(i-1),x(i-1)+h));