Matlab-10:Ritz-Galerkin方法求解二阶常微分方程

一、代数多项式法：

 1 tic;
 2 clear
 3 clc
 4 % N=input('please key in the value of ''N''');
 5 N=10;
 6 M=100;
 7 h=1/M;
 8 X=0:h:1;
 9 accurate_fun=inline('x.^2 - (2*exp(x))/(exp(1) + 1) - (2*exp(-x)*exp(1))/(exp(1) + 1) + 2');
10  f=inline('x.^2-x');
11  phi=inline('x.*(1-x).*x.^(i-1)','i','x');
12  diff_phi=inline('i*x.^(i-1)-(i+1)*x.^i','i','x');
13  for j=1:N
14      for i=1:N
15 A(i,j)=quad(@(x)phi(i,x).*phi(j,x)+diff_phi(i,x).*diff_phi(j,x),0,1);
16      end
17      b(j,1)=quad(@(x) phi(j,x).*f(x),0,1);
18  end
19  C=A;
20 syms x;
21  Un=0;
22 for i=1:N
23 Un=Un+C(i)*phi(i,x);
24 end
25 Un=Un+x;
26  numerical= double(vpa(subs(Un,'x',X)));
27  accurate=accurate_fun(X);
28  subplot(1,2,1)
29  plot(X,numerical,'r -',X,accurate,'b >');
30   title('numerical VS accurate');
31  legend('numerical','accurate');
32  grid on;
33   subplot(1,2,2);
34   plot(X,numerical-accurate,'g');
35   title('error');
36   grid on;
37  toc;

二、三角函数法：

 1 tic;
 2 clear
 3 clc
 4 % N=input('please key in the value of ''N''');
 5 N=10;
 6 M=100;
 7 h=1/M;
 8 X=0:h:1;
 9 accurate_fun=inline('x.^2 - (2*exp(x))/(exp(1) + 1) - (2*exp(-x)*exp(1))/(exp(1) + 1) + 2');
10  f=inline('x.^2-x');
11  phi=inline('sin(i*pi*x)','i','x');
12  diff_phi=inline('i*pi*cos(i*pi*x)','i','x');
13  for j=1:N
14      for i=1:N
15 A(i,j)=quad(@(x)phi(i,x).*phi(j,x)+diff_phi(i,x).*diff_phi(j,x),0,1);
16      end
17      b(j,1)=quad(@(x) phi(j,x).*f(x),0,1);
18  end
19  C=A;
20  syms x;
21  Wn=0;
22 for i=1:N
23 Wn=Wn+C(i)*phi(i,x);
24 end
25 Un=Wn+x;
26 numerical=vpa(subs(Un,'x',X));
27 accurate=accurate_fun(X);
28 subplot(1,2,1)
29 plot(X,numerical,'r -',X,accurate,'g ^');
30 title('Ritz Galerkin method');
31 legend('numerical solution','accurate solution');
32 grid on;
33 subplot(1,2,2)
34 plot(X,numerical-accurate,'b');
35 title('error');
36 grid on;
37 toc;