Matlab-7:偏微分方程数值解法-李荣华-有限元解导数边界值的常微分（Galerkin方法）

p47.（实习题-李荣华）用线性元求下列边值问题的数值解

 tic;

 % this method is transform from Galerkin method

 %also call it as finit method

 %is used for solving two point BVP which is the first and second term.

 %this code was writen by HU.D.dong in February 11th

 %MATLAB 7.0

 clear;

 clc;

 N=;

 h=/N;

 X=:h:;

 f=inline('(0.5*pi^2)*sin(0.5*pi.*x)');

 %以下是右端向量:

 for i=:N

         fun1=@(x) pi^/.*sin(pi/.*x).*(-(x-X(i))/h);

          fun2=@(x) pi^/.*sin(pi/.*x).*((x-X(i-))/h);

     f_phi(i-,)=quad(fun1,X(i),X(i+))+quad(fun2,X(i-),X(i));

     end

  funN=@(x) pi^/.*sin(pi/.*x).*(x-X(N))/h;

     f_phi(N)=quad(funN,X(N),X(N+));

     %以下是刚度矩阵：

  A11=quad(@(x) /h+0.25*pi^*h.*(-*x+*x.^),,);

  A12=quad(@(x) -/h+0.25*pi^*h.*(-x).*x,,);

  ANN=quad(@(x) /h+0.25*pi^*h*x.^,,);

  A=diag([A11*ones(,N-),ANN],)+diag(A12*ones(,N-),)+diag(A12*ones(,N-),-);

  Numerical_solution=A\f_phi;

  Numerical_solution=[;Numerical_solution];

     %Accurate solution on above以下是精确解

     %%

     for i=:length(X)

    Accurate_solution(i,)=sin((pi*X(i))/)/ - cos((pi*X(i))/)/ + exp((pi*X(i))/)*((exp(-(pi*X(i))/)*cos((pi*X(i))/))/ + (exp(-(pi*X(i))/)*sin((pi*X(i))/))/);

     end

     figure();

     grid on;

     subplot(,,);

      plot(X,Numerical_solution,'ro-',X,Accurate_solution,'b^:');

     title('Numerical solutions vs Accurate solutions');

     legend('Numerical_solution','Accurate_solution');

     subplot(,,);

       plot(X,Numerical_solution-Accurate_solution,'b x');

     legend('error_solution');

     title('error');

     toc;