Matlab:导数边界值的有限元(Ritz)法
1 tic; 2 % this method is transform from Ritz method 3 %is used for solving two point BVP 4 %this code was writen by HU.D.dong in February 11th 2017 5 %MATLAB 7.0 6 clear 7 clc 8 N=50; 9 h=1/N; 10 X=0:h:1; 11 f=inline('pi^2/2*sin(pi/2*x)'); 12 %以下是右端向量； 13 for i=2:N 14 fun1=@(x) f(X(i-1)+h*x).*x+f(X(i)+h*x).*(1-x); 15 fi(i-1,1)=h*quad(fun1,0,1); 16 end 17 funN=@(x) f(X(N-1)+h*x).*x; 18 fi(N,1)=h*quad(funN,0,1); 19 %以下是stiff矩阵； 20 a11=1/h+pi^2*h/12; 21 aii=2*a11; 22 aij=-1/h+pi^2*h/24; 23 A=diag(aii*ones(N,1),0)+diag(aij*ones(N-1,1),1)+diag(aij*ones(N-1,1),-1); 24 A(N,N)=a11; 25 numerical_solution=A\fi; 26 numerical_solution=[0;numerical_solution]; 27 %以下是真解； 28 for i=1:length(X) 29 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); 30 end 31 grid on; 32 subplot(1,2,1); 33 plot(X,numerical_solution,'ro-',X,Accurate_solution,'b^:'); 34 title('Numerical solutions vs Accurate solutions'); 35 legend('Numerical_solution','Accurate_solution'); 36 subplot(1,2,2); 37 plot(X,numerical_solution-Accurate_solution,'b x'); 38 legend('error_solution'); 39 title('error'); 40 toc;
效果图：
胡冬冬
