110,566
社区成员
发帖
与我相关
我的任务
分享高斯消元 求四元！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！
高斯消元 求四元
na0 + a1∑xi + a2∑xi^2 + a3∑xi^3 = ∑yi
a0∑xi + a1∑xi^2 + a2∑xi^3 + a3∑xi^4 = ∑yixi
a0∑xi^2 + a1∑xi^3 + a2∑xi^4 + a3∑xi^5 = ∑yixi^2
a0∑xi^3 + a1∑xi^4 + a2∑xi^5 + a3∑xi^6 = ∑yixi^3
x，y为坐标，已知
n为点数，已知
i已知
求： a0,a1,a2,a3
以下是目前的做法，自己也不知道对不对。。。。。。。根据人家的三元改的。。请帮看下。
na0 + a1∑xi + a2∑xi^2 + a3∑xi^3 = ∑yi
a0∑xi + a1∑xi^2 + a2∑xi^3 + a3∑xi^4 = ∑yixi
a0∑xi^2 + a1∑xi^3 + a2∑xi^4 + a3∑xi^5 = ∑yixi^2
a0∑xi^3 + a1∑xi^4 + a2∑xi^5 + a3∑xi^6 = ∑yixi^3
x，y为坐标，已知
n为点数，已知
i已知
求： a0,a1,a2,a3
以下是目前的做法，自己也不知道对不对。。。。。。。根据人家的三元改的。。请帮看下。
//0：na0 + a1∑xi + a2∑xi^2 + a3∑xi^3 = ∑yi
//1: a0∑xi + a1∑xi^2 + a2∑xi^3 + a3∑xi^4 = ∑yixi
//2：a0∑xi^2 + a1∑xi^3 + a2∑xi^4 + a3∑xi^5 = ∑yixi^2
//3：a0∑xi^3 + a1∑xi^4 + a2∑xi^5 + a3∑xi^6 = ∑yixi^3
public static void Main(string[] args) // 主函数
{
int yuan = 4;
int n = 2;
DataTable dt = new DataTable("PointF");
dt.Columns.Add("x");
dt.Columns.Add("y");
DataRow dr = dt.NewRow();
dr["x"] = 460;
dr["y"] = 0.1107;
dt.Rows.Add(dr);
dr = dt.NewRow();
dr["x"] = 480;
dr["y"] = 0.1146;
dt.Rows.Add(dr);
double[] xx = new double[dt.Rows.Count];
double[] y = new double[dt.Rows.Count];
double[] xy1 = new double[dt.Rows.Count];
double[] xy2 = new double[dt.Rows.Count];
double[] xy3 = new double[dt.Rows.Count];
for (int i = 0; i < dt.Rows.Count; i++)
{
xx[i] = double.Parse(dt.Rows[i]["x"].ToString());
y[i] = double.Parse(dt.Rows[i]["y"].ToString());
xy1[i] = double.Parse(dt.Rows[i]["y"].ToString()) * double.Parse(dt.Rows[i]["y"].ToString());
xy2[i] = Math.Pow(double.Parse(dt.Rows[i]["y"].ToString()), 2) * double.Parse(dt.Rows[i]["y"].ToString());
xy3[i] = Math.Pow(double.Parse(dt.Rows[i]["y"].ToString()), 3) * double.Parse(dt.Rows[i]["y"].ToString());
}
//转换矩阵
double[,] jz =
{
{ n,GetSGM(xx, 1), Math.Pow(GetSGM(xx, 1), 2), Math.Pow(GetSGM(xx, 1), 3),GetSGM(y, 1)}
, { GetSGM(xx, 1), Math.Pow(GetSGM(xx, 1), 2), Math.Pow(GetSGM(xx, 1), 3), Math.Pow(GetSGM(xx, 1), 4), GetSGM(xy1, 1) }
, { Math.Pow(GetSGM(xx, 1),2), Math.Pow(GetSGM(xx, 1),3), Math.Pow(GetSGM(xx, 1),4), Math.Pow(GetSGM(xx, 1),5),GetSGM(xy2, 1) }
, { Math.Pow(GetSGM(xx, 1),3), Math.Pow(GetSGM(xx, 1),4), Math.Pow(GetSGM(xx, 1),5), Math.Pow(GetSGM(xx, 1),6),GetSGM(xy3, 1) }
};
double[] a = new double[yuan];
Gauss(yuan, jz, a);
// 输出方程组的解
Console.WriteLine("方程组的解为：");
for (int i = 0; i < yuan; i++)
{
Console.Write("x({0})={1} ", i, a[i]);
}
Console.WriteLine();
}
public static void Gauss(int n, double[,] a, double[] x)
{
double d;
Console.WriteLine("高斯消去法解方程组的中间过程");
Console.WriteLine("============================");
Console.WriteLine("中间过程");
Console.WriteLine("增广矩阵：");
printArray(n, a); Console.WriteLine();
// 消元
for (int k = 0; k < n; k++)
{
Console.WriteLine("第{0}步", k + 1);
Console.WriteLine("初始矩阵：");
printArray(n, a); Console.WriteLine();
selectMainElement(n, k, a); // 选择主元素
Console.WriteLine("选择主元素后的矩阵：");
printArray(n, a); Console.WriteLine();
d = a[k, k];
for (int j = k; j <= n; j++)
{
a[k, j] = a[k, j] / d;
}
Console.WriteLine("将第{0}行中a[{0},{0}]化为1后的矩阵：", k + 1);
printArray(n, a); Console.WriteLine();
for (int i = k + 1; i < n; i++)
{
d = a[i, k]; // 这里使用变量d将a[i,k]的值保存下来的原理与上面注释中说明的一样
for (int j = k; j <= n; j++)
{
a[i, j] = a[i, j] - d * a[k, j];
}
}
Console.WriteLine("消元后的矩阵：");
printArray(n, a);
Console.WriteLine();
x[n - 1] = a[n - 1, n];
for (int i = n - 1; i >= 0; i--)
{
x[i] = a[i, n];
for (int j = i + 1; j < n; j++) {
x[i] = x[i] - a[i, j] * x[j];
}
}
}
}
public static void selectMainElement(int n, int k, double[,] a)
{
double t, mainElement; // mainElement用于保存主元素的值
int l; // 用于保存主元素所在的行号
mainElement = Math.Abs(a[k, k]); // 注意别忘了取绝对值
l = k;
for (int i = k + 1; i < n; i++)
{
if (mainElement < Math.Abs(a[i, k]))
{
mainElement = Math.Abs(a[i, k]);
l = i; // 记下主元素所在的行号
}
}
// l是主元素所在的行。将l行与k行交换，每行前面的k个元素都是0，不必交换
if (l != k)
{
for (int j = k; j <= n; j++)
{
t = a[k, j]; a[k, j] = a[l, j]; a[l, j] = t;
}
}
}
public static void printArray(int n, double[,] a)
{
for (int i = 0; i < n; i++)
{
for (int j = 0; j <= n; j++)
{
Console.Write("{0,10:F8} ", a[i, j]);
//Console.Write(a[i, j]);
}
Console.WriteLine();
}
}
//求 和：∑
public static double GetSGM(double[] x,int count)
{
double a = 0;
for (int i = 0; i < count; i++)
{
a += x[i];
}
//此处不除1000就溢出了。
return a/1000;
}
...全文
请发表友善的回复…
发表回复
gwf25sz 2012-02-07
- 打赏
- 举报
继续求助
gwf25sz 2012-02-03
- 打赏
- 举报
没了？？？？？？？？？？？
gwf25sz 2012-02-02
- 打赏
- 举报
[Quote=引用 3 楼 hhqsy 的回复:]
换更高精度的数据类型看看
[/Quote]
精度倒不是问题，我想知道，我现在的公式，用我的代码做，方法对不对。。。。。
换更高精度的数据类型看看
[/Quote]
精度倒不是问题，我想知道，我现在的公式，用我的代码做，方法对不对。。。。。
gomoku 2012-02-02
- 打赏
- 举报
高斯消元的原理很直接（用纸笔很容易做），不过要写一个好程序很不容易。
因为数值计算本身有误差，如何控制误差的传递和扩大，并不是很容易做好的。
"Numerical Recipes"是一本非常好的书，也有源码。如果有这方面的爱好，很值得参考。
如果不是兴趣所在，下载一些免费库是更好的选择。
因为数值计算本身有误差，如何控制误差的传递和扩大，并不是很容易做好的。
"Numerical Recipes"是一本非常好的书，也有源码。如果有这方面的爱好，很值得参考。
如果不是兴趣所在，下载一些免费库是更好的选择。
昵称不能为空噢噢噢噢 2012-02-02
- 打赏
- 举报
换更高精度的数据类型看看
昵称不能为空噢噢噢噢 2012-02-02
- 打赏
- 举报
decimal
