首页 > 代码库 > Matlab:Ritz-Galerkin方法求解二阶常微分方程
Matlab: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\b; 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\b; 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;
Matlab:Ritz-Galerkin方法求解二阶常微分方程
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。