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;