cicyabc
2007-06-14, 12:18
本人毕业设计涉及了数值模拟,于是自学了一点matlab用于求解线性方程组,导入矩阵各元素后发现矩阵的秩只有61,换了好多次边界条件,秩仍为61虽能求解,但是结果明显有误。不知道是不是程序有问题。 前面的程序很短的,望好心人抽空给于指正,万分感谢。
下面是前面的一部分程序:
%根据现有资料,计算下一时刻的水位,流量,盐度
A0=zeros(62,62);
format long e ;
T1=load('t1.txt');%T1=Δt/bi/Δx
T2=load('t2.txt');%T2=ΔtgA/Δx
T3=load('t3.txt');%T3=1-ΔtUi/Δx
T4=load('t4.txt');%T4=ΔtU(i+1)/Δx
T5=load('t5.txt');%T5=-ΔtgAdc/Δx/(Si+0.741)
T6=load('t6.txt');%T6=1-αΔt/Δx +2KxΔt/Δx
T7=load('t7.txt');%T7=αΔt/Δx-KxΔt/Δx/Δx
T8=load('t8.txt');%T8=KxΔt/Δx/Δx
%i=1时方程(1)的系数
A0(1,1)=1;
A0(1,22)=T1(1);
%i=21时方程(1)的系数
A0(21,21)=1;
A0(21,41)=-T1(21);
A0(21,42)=T1(21);
%i=1时方程(2)的系数
A0(22,1)=-T2(1);
A0(22,2)=T2(1);
A0(22,22)=T4(1);
A0(22,43)=-T5(1);
%i=21时方程(2)的系数
A0(42,21)=-T2(21);
A0(42,41)=T3(21);
A0(42,42)=T4(21);
A0(42,62)=T5(21);
%i=2时方程(3)的系数
A0(43,43)=T6(2);
A0(43,44)=T7(2);
%i=21时方程(3)的系数
A0(62,61)=T8(21);
A0(62,62)=T6(21);
%根据方程(1)(2)对其他系数赋值
for i=2:20;
A0(i,i)=1;
A0(i,i+20)=-T1(i);
A0(i,i+21)=T1(i);
A0(i+21,i)=-T2(i);
A0(i+21,i+1)=T2(i);
A0(i+21,i+20)=T3(i);
A0(i+21,i+21)=T4(i);
A0(i+21,i+41)=T5(i);
A0(i+21,i+42)=-T5(i);
end;
%根据方程(3)对其他系数赋值
for i=3:20;
A0(i+41,i+41)=T6(i);
A0(i+41,i+42)=T7(i);
A0(i+41,i+40)=T8(i);
end;
B0=load('方程右边.txt');
X1=A0\B0;%第一次求解结束
X1,%输出X1
运行是总出现一下警告:
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.721593e-023.
下面是前面的一部分程序:
%根据现有资料,计算下一时刻的水位,流量,盐度
A0=zeros(62,62);
format long e ;
T1=load('t1.txt');%T1=Δt/bi/Δx
T2=load('t2.txt');%T2=ΔtgA/Δx
T3=load('t3.txt');%T3=1-ΔtUi/Δx
T4=load('t4.txt');%T4=ΔtU(i+1)/Δx
T5=load('t5.txt');%T5=-ΔtgAdc/Δx/(Si+0.741)
T6=load('t6.txt');%T6=1-αΔt/Δx +2KxΔt/Δx
T7=load('t7.txt');%T7=αΔt/Δx-KxΔt/Δx/Δx
T8=load('t8.txt');%T8=KxΔt/Δx/Δx
%i=1时方程(1)的系数
A0(1,1)=1;
A0(1,22)=T1(1);
%i=21时方程(1)的系数
A0(21,21)=1;
A0(21,41)=-T1(21);
A0(21,42)=T1(21);
%i=1时方程(2)的系数
A0(22,1)=-T2(1);
A0(22,2)=T2(1);
A0(22,22)=T4(1);
A0(22,43)=-T5(1);
%i=21时方程(2)的系数
A0(42,21)=-T2(21);
A0(42,41)=T3(21);
A0(42,42)=T4(21);
A0(42,62)=T5(21);
%i=2时方程(3)的系数
A0(43,43)=T6(2);
A0(43,44)=T7(2);
%i=21时方程(3)的系数
A0(62,61)=T8(21);
A0(62,62)=T6(21);
%根据方程(1)(2)对其他系数赋值
for i=2:20;
A0(i,i)=1;
A0(i,i+20)=-T1(i);
A0(i,i+21)=T1(i);
A0(i+21,i)=-T2(i);
A0(i+21,i+1)=T2(i);
A0(i+21,i+20)=T3(i);
A0(i+21,i+21)=T4(i);
A0(i+21,i+41)=T5(i);
A0(i+21,i+42)=-T5(i);
end;
%根据方程(3)对其他系数赋值
for i=3:20;
A0(i+41,i+41)=T6(i);
A0(i+41,i+42)=T7(i);
A0(i+41,i+40)=T8(i);
end;
B0=load('方程右边.txt');
X1=A0\B0;%第一次求解结束
X1,%输出X1
运行是总出现一下警告:
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.721593e-023.