左除怎么求不出正确的线性代数方程组的解！重点部分标示

[meerbao]

左除怎么求不出正确的线性代数方程组的解！重点部分标示
%%-----------------
%-说明，本程序所有的循环应该没有问题，
%-主要是利用laplace的延迟性质，求解地层的压力导数与时间的关系
%-主要遇到的问题是在，计算F时得到的F不是理想的结果
%-请各位高手帮下忙，帮我看看到底是什么原因
%-pp.txt和tt.txt文件在附件里，如果必要麻烦高手调试下，谢谢！！
%%-----------------

clear;
clc;
%function [a b c]=Delaytime2()
p=textread('pp.txt',' %12.4f\r\n');%压差
tt=textread('tt.txt',' %12.4f\r\n');%时间
%M=length(tt)-20;

PP = p(1:12);%%做了修改，防止矩阵g奇异
TT = tt(1:12);
PD=1.0*30.48*PP/(1.842*10^(-3)*1462.683*1.24*1.24); %仅仅是变量的转换PP->PD
tD=3.6*1.0*TT/0.27/0.0015954/1.24/0.01/24; %仅仅是变量的转换TT->tD
M = length(PP);
NS = M; %必须是偶数

Vt = zeros(NS,1);
NN=NS/2;
for i = 1:NS
VV=0;
k1 = floor((i+1)/2);
k2 = min(i,floor(NN));
for k = k1:k2
Fz = k^(NN)*factorial(2*k);
Fm = factorial(k-1)*factorial(k)*factorial((NN)-k)*factorial(i-k)*factorial(2*k-i);
VV = VV + Fz/Fm;
end
VV=(-1)^(NN+i)*VV;
Vt(i)=VV;
end
dlmwrite('m1.txt',Vt,'delimiter','\t','precision', 20);%输出Vt
V=dlmread('m1.txt');
%% Vt-V %%已验证结果为0
g = zeros(M,M);
for i = 1:M
for j =1: M
g(i,j) = V(j)*exp(-log(2)*(1-tD(i)/tD(M))*j);%%得到M*M维矩阵
end
end

dlmwrite('m2.txt',g,'delimiter','\t','precision', 6);%输出g
G=dlmread('m2.txt');
%G-g %%已验证结果为0
%det(G) %已验证结果不是0，也不是inf

R = PD*tD(M)/log(2);
dlmwrite('m5.txt',R,'delimiter','\t','precision', 6);%输出R
r=dlmread('m5.txt');
%R-r %%已验证结果为0

F = G\r;%左除计算F

% G*F-r %%结果不是0
%------不知道为什么会这样？？？？？？
%[L,U]=lu(G);%%L-U分解计算F
%F=U\(L\r);
dlmwrite('m3.txt',F,'delimiter','\t','precision', 6);% 输出F
L=dlmread('m3.txt');

for i=1:M
ff=zeros(M,1);
fV=0;
for j=1:M
ff(j)=ff(j)+V(j)*exp(-log(2)*(1-tD(i)/tD(M))*j)*(log(2)/tD(M)*j)*L(j);
end
for j=1:M
fV=fV+ff(j);
end
f(i)=fV;
end
f=f';
dlmwrite('m4.txt',f,'delimiter','\t','precision', 6);
ef=dlmread('m4.txt');
df=ef*log(2)/tD(M); %%压力导数
%df=ef*log(2); %%压力导数*时间

figure();
plot(tD,df,'ro');
xlabel('tD');
ylabel('pD*tD');
title('压差导数乘时间--时间');

[slgu]

试试下列小程序：
x=[0.0032, 0.0074, 0.0107, 0.0136, 0.0166, 0.0224, 0.0282, 0.0341, 0.0399, 0.0457, 0.0545, 0.0632, 0.072, 0.0807, 0.0903, 0.1102, 0.131, 0.1557, 0.1807, 0.2107, 0.2449, 0.2852, 0.3391, 0.4204, 0.5282, 0.591, 0.671, 0.7899, 0.9693, 1.086, 1.2282, 1.4004, 1.6016, 1.8382, 2.1216, 2.466, 2.8771, 3.3882, 3.9971, 4.7438, 5.6816, 6.8193, 7.9749]';
y=[0.04, 0.07, 0.1, 0.13, 0.16, 0.22, 0.28, 0.34, 0.41, 0.47, 0.56, 0.65, 0.72, 0.82, 0.91, 1.09, 1.27, 1.46, 1.63, 1.82, 2.0, 2.17, 2.35, 2.53, 2.67, 2.73, 2.79, 2.85, 2.91, 2.94, 2.97, 3.0, 3.03, 3.06, 3.09, 3.12, 3.15, 3.18, 3.21, 3.21, 3.27, 3.3, 3.33]';
fx=@(b,x)(b(1)*x+b(2)*x.^2+b(3)*x.^3)./(1+b(4)*x+b(5)*x.^2+b(6)*x.^3);
b=[9.7865318, 32.775392, 27.227736, 1.967689, 13.118039, 7.8103828];
x1=linspace(0,8)
plot(x,y,'o','markerfacecolor','k','linewidth',3);
hold on
y1=fx(b,x1);
plot(x1,y1,'linewidth',2);