【发布时间】：2016-09-27 05:32:40
【问题描述】：
我正在尝试使用 Thomas 算法求解微分热方程。
物理问题：我们有插头，左侧温度为0，右侧温度为1。
对于 Thomas 算法，我编写了一个函数，它接受三个 QVector 和 int 值数量的方程。
这是我的代码：
#include <QCoreApplication>
#include <QVector>
#include <QDebug>
#include <iostream>
using std::cin;
void enterIn(QVector<float> &Array, int Amount_of_elements){
int transit;
for(int i=0;i<Amount_of_elements;i++){
cin>>transit;
Array.push_back(transit);
}
}
QVector<float> shuttle_method(const QVector<float> &below_main_diagonal,
QVector<float> &main_diagonal,
const QVector<float> &beyond_main_diagonal,
const QVector<float> &free_term,
const int N){
QVector <float> c;
QVector <float> d;
for(int i=0;i<N;i++){
main_diagonal[i]*=(-1);
}
QVector<float> x; //result
c.push_back(beyond_main_diagonal[0]/main_diagonal[0]);
d.push_back(-free_term[0]/main_diagonal[0]);
for(int i=1;i<=N-2;i++){
c.push_back(beyond_main_diagonal[i]/(main_diagonal[i]-below_main_diagonal[i]*c[i-1]));
d.push_back( (below_main_diagonal[i]*d[i-1] - free_term[i]) / (main_diagonal[i]- below_main_diagonal[i]*c[i-1]) );
}
x.resize(N);
//qDebug()<<x.size()<<endl;
int n=N-1;
x[n]=(below_main_diagonal[n]*d[n-1]-free_term[n])/(main_diagonal[n]-below_main_diagonal[n]*c[n-1]);
for(int i=n-1;i>=0;i--){
x[i]=c[i]*x[i+1]+d[i];
// qDebug()<<x[i]<<endl;
}
return x;
}
int main()
{
QVector <float> alpha; // below
QVector <float> beta; // main diagonal * (-1)
QVector <float> gamma; // beyond
QVector <float> b; // free term
QVector<float> T;
int cells_x=40; //amount of equations
alpha.resize(cells_x);
beta.resize(cells_x);
gamma.resize(cells_x);
b.resize(cells_x);
T.resize(cells_x);
float dt=0.2,h=0.1;
alpha[0]=0;
for(int i=1;i<cells_x;i++){
alpha[i]= -dt/(h*h);
}
for(int i=0;i<cells_x;i++){
beta[i] = (2*dt)/(h*h)+1;
}
for(int i=0;i<cells_x-1;i++){
gamma[i]= -dt/(h*h);
}
gamma[cells_x-1]=0;
qDebug()<<"alpha= "<<endl<<alpha.size()<<alpha<<endl<<"beta = "<<endl<<beta.size()<<beta<<endl<<"gamma= "<<gamma.size()<<gamma<<endl;
for(int i=0;i<cells_x-1;i++){
T[i]=0;
}
T[cells_x-1]=1;
qDebug()<<endl<<endl<<T<<endl;
//qDebug()<< shuttle_method(alpha,beta,gamma,b,N);
QVector<float> Tn;
Tn.resize(cells_x);
Tn = shuttle_method(alpha,beta,gamma,T,cells_x);
Tn[0]=0;Tn[cells_x-1]=1;
for(int stepTime = 0; stepTime < 50; stepTime++){
Tn = shuttle_method(alpha,beta,gamma,Tn,cells_x);
Tn[0]=0;
Tn[cells_x-1]=1;
qDebug()<<Tn<<endl;
}
return 0;
}
我的问题是： 当我编译并运行它时，我得到了这个：
Tn <20 items> QVector<float>
0 float
0.000425464 float
0.000664658 float
0.000937085 float
0.00125637 float
0.00163846 float
0.00210249 float
0.00267163 float
0.00337436 float
0.00424581 float
0.00532955 float
0.00667976 float
0.00836396 float
0.0104664 float
0.0130921 float
0.0163724 float
0.0204714 float
0.0255939 float
0.0319961 float
Tn <20 items> QVector<float>
0 float
-0.000425464 float
0.000643385 float
-0.000926707 float
0.00120951 float
-0.00161561 float
0.00202056 float
-0.00263167 float
0.00324078 float
-0.00418065 float
0.00511726 float
-0.00657621 float
0.00802998 float
-0.0103034 float
0.0125688 float
-0.0161171 float
0.0196527 float
-0.0251945 float
0.0307164 float
1 float
Tn <20 items> QVector<float>
0 float
0.000425464 float
0.000664658 float
0.000937085 float
0.00125637 float
0.00163846 float
0.00210249 float
0.00267163 float
0.00337436 float
0.00424581 float
0.00532955 float
0.00667976 float
0.00836396 float
0.0104664 float
0.0130921 float
0.0163724 float
0.0204714 float
0.0255939 float
0.0319961 float
Tn <20 items> QVector<float>
0 float
-0.000425464 float
0.000643385 float
-0.000926707 float
0.00120951 float
-0.00161561 float
0.00202056 float
-0.00263167 float
0.00324078 float
-0.00418065 float
0.00511726 float
-0.00657621 float
0.00802998 float
-0.0103034 float
0.0125688 float
-0.0161171 float
0.0196527 float
-0.0251945 float
0.0307164 float
1 float
一次又一次循环。
我不知道为什么我会得到这个。
也许我的错误在于分配 Tn 我的 Thomas-method-function 的结果？
或实现托马斯方法？还是在边界条件下？
【问题讨论】：
-
您是否尝试过在调试器中运行和单步执行您的代码？这通常会有所帮助。从小数字开始，这样循环就可以快速完成。
-
调试器是解决此类问题的正确工具。 在询问 Stack Overflow 之前，您应该逐行浏览您的代码。如需更多帮助，请阅读How to debug small programs (by Eric Lippert)。至少，您应该 [编辑] 您的问题，以包含一个重现您的问题的 Minimal, Complete, and Verifiable 示例，以及您在调试器中所做的观察。
-
我当然用过！但这没有帮助，我仍然得到相同的结果，但不知道为什么
-
为什么要否定主对角线作为函数的第一步？然后创建一个三对角矩阵类并使用双精度而不是浮点数会很有用。