P47.（实习题-李荣华）用线性元求下列边值问题的数值解
1 tic; 2 % this method is transform from Galerkin method 3 %also call it as finit method 4 %is used for solving two point BVP which is the first and second term. 5 %this code was writen by HU.D.dong in February 11th 2017 6 %MATLAB 7.0 7 clear; 8 clc; 9 N=50; 10 h=1/N; 11 X=0:h:1; 12 f=inline(\'(0.5*pi^2)*sin(0.5*pi.*x)\'); 13 %以下是右端向量: 14 for i=2:N 15 fun1=@(x) pi^2/2.*sin(pi/2.*x).*(1-(x-X(i))/h); 16 fun2=@(x) pi^2/2.*sin(pi/2.*x).*((x-X(i-1))/h); 17 f_phi(i-1,1)=quad(fun1,X(i),X(i+1))+quad(fun2,X(i-1),X(i)); 18 end 19 funN=@(x) pi^2/2.*sin(pi/2.*x).*(x-X(N))/h; 20 f_phi(N)=quad(funN,X(N),X(N+1)); 21 %以下是刚度矩阵： 22 A11=quad(@(x) 2/h+0.25*pi^2*h.*(1-2*x+2*x.^2),0,1); 23 A12=quad(@(x) -1/h+0.25*pi^2*h.*(1-x).*x,0,1); 24 ANN=quad(@(x) 1/h+0.25*pi^2*h*x.^2,0,1); 25 A=diag([A11*ones(1,N-1),ANN],0)+diag(A12*ones(1,N-1),1)+diag(A12*ones(1,N-1),-1); 26 Numerical_solution=A\f_phi; 27 Numerical_solution=[0;Numerical_solution]; 28 %Accurate solution on above以下是精确解 29 %% 30 for i=1:length(X) 31 Accurate_solution(i,1)=sin((pi*X(i))/2)/2 - cos((pi*X(i))/2)/2 + exp((pi*X(i))/2)*((exp(-(pi*X(i))/2)*cos((pi*X(i))/2))/2 + (exp(-(pi*X(i))/2)*sin((pi*X(i))/2))/2); 32 end 33 figure(1); 34 grid on; 35 subplot(1,2,1); 36 plot(X,Numerical_solution,\'ro-\',X,Accurate_solution,\'b^:\'); 37 title(\'Numerical solutions vs Accurate solutions\'); 38 legend(\'Numerical_solution\',\'Accurate_solution\'); 39 subplot(1,2,2); 40 plot(X,Numerical_solution-Accurate_solution,\'b x\'); 41 legend(\'error_solution\'); 42 title(\'error\'); 43 toc;
.（实习题）用线性元求下列边值问题的数值解
解：从出发可知：
形式的变分方程为
,
其中
,
又
因此
在单元中,应用仿射变换(局部坐标)
节点基函数为
由此可知系数矩阵是一个三对角矩阵且对称的矩阵：