Labfans是一个针对大学生、工程师和科研工作者的技术社区。 | 论坛首页 | 联系我们(Contact Us) |
![]() |
![]() |
#1 |
初级会员
注册日期: 2009-06-08
帖子: 1
声望力: 0 ![]() |
![]()
左除怎么求不出正确的线性代数方程组的解!重点部分标示
%%----------------- %-说明,本程序所有的循环应该没有问题, %-主要是利用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('压差导数乘时间--时间'); 此帖于 2009-06-08 10:29 被 meerbao 编辑。 |
![]() |
![]() |