问题原型
边界条件:
(1)取 alpha,beta=1 f=x
(2)p r q=1
(4)(5)
实际结果如下:
完整代码如下
clc; clear all
M=40; % Local-segment number
N=M+1;%points number
av=1;bv=1;L=1;le=L/M;
alpha=ones(M,1)'*av;
beta=ones(M,1)'*bv;
l=ones(M,1)'*le;
f=[1:M]*le;
% 2. boundary condition
p=1;%x=0 position
r=1;%x=L position
q=1;
%3 print input
%4 get K: K(i,i) to a , k(i+1,i) to c
i=2:N-1;
a=[alpha(1)/l(1)+beta(1)*l(1)/3,...
alpha(i-1)./l(i-1)+beta(i-1).*l(i-1)/3 ...
+ alpha(i)./l(i)+beta(i).*l(i)/3, ...
alpha(M)/l(M)+beta(M)*l(M)/3]';
i=1:N-1;
c=[-alpha(i)./l(i)+beta(i).*l(i)/6]';
%5 get b
i=2:N-1;
b=[f(1)*l(1)/2, ...
f(i-1).*l(i-1)/2+f(i).*l(i)/2, ...
f(M)*l(M)/2]';
%6 boudary condition apply
a(end)=a(end)+r; %x=L third condition
b(end)=b(end)+q;
K=diag(a,0)+diag(c,-1)+diag(c,+1);
K(1,1)=1; b(1)=p;K(1,2:end)=0;%x=0 1st condition
%7 compute
x=(0:M)*le;
y=K^(-1)*b;
plot(x,y);
hold on;
x1=(0:100)*0.01;
c1=-1/2/exp(1);
c2=1-c1;
y1=c1*exp(x1)+c2*exp(-x1)+x1;
plot(x1,y1,'-r');
hold off;