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=Af_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;
效果图：
