一、最小二乘法,用MATLAB实现
1. 数值实例
下面给定的是乌鲁木齐最近1个月早晨7:00左右(新疆时间)的天气预报所得到的温度,按照数据出任意次曲线拟合方程和它的图像。下面用MATLAB编程对上述数据进行最小二乘拟合。
下面用MATLAB编程对上述数据进行最小二乘拟合
2、程序代码
x=[1:1:30];
y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1];
a1=polyfit(x,y,3) %三次多项式拟合%
a2= polyfit(x,y,9) %九次多项式拟合%
a3= polyfit(x,y,15) %十五次多项式拟合%
b1=polyval(a1,x)
b2=polyval(a2,x)
b3=polyval(a3,x)
r1= sum((y-b1).^2) %三次多项式误差平方和%
r2= sum((y-b2).^2) %九次次多项式误差平方和%
r3= sum((y-b3).^2) %十五次多项式误差平方和%
plot(x,y,'*') %用*画出x,y图像%
hold on
plot(x,b1, 'r') %用红线画出x,b1图像%
hold on
plot(x,b2, 'g') %用绿线画出x,b2图像%
hold on
plot(x,b3, 'b:o') %用蓝o线画出x,b3图像%
3、数值结果
不同次数多项式拟合误差平方和为:
r1=67.6659
r2=20.1060
r3=3.7952
r1、r2、r3分别表示三次、九次、十五次多项式误差平方和。
4、拟合曲线如下图
二、 线性方程组的求解( 高斯-塞德尔迭代算法 )
1、实例: 求解线性方程组(见书P233页)
⎪⎪⎩⎪⎪⎨⎧=++=-+=+-36
12363311420238321321321x x x x x x x x x  记A x=b, 其中
⎥⎥
⎥⎦
⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦
⎢⎢⎢⎣⎡=⎥⎥
⎥⎦
⎢⎢⎢⎣⎡--=363320,
,12361114238321b x A x x x
任取初始值
()
()T
x
0000=,进行迭代。
2、用C 语言实现,源程序如下:
#include<stdio.h>
#include<stdlib.h>
struct Node
{
double value;  // 矩阵非零元素值
int index;  // 元素的索引(序号或列标)
};
struct Matrix
{
struct Node *data; // 记录元素数据
int total_elem; // 矩阵非零元素个数
int total_ln; // 矩阵的总行数
int total_col; // 矩阵的总列数
};
// 采用MSR法实现系数矩阵的存储
void Create_Matrix( struct Matrix &matrix )
{
int index = 0, i = 0;
double value = 0.0;
printf( "请输入线性方程组系数矩阵的非零元素个数及总行数、总列数: " ); scanf( "%d%d%d", &al_elem, &al_ln, &al_col );
matrix.data = (struct Node *)malloc( sizeof(struct Node)*(al_elem+1) );
printf( "请依次输入矩阵的主对角线元素及每行的非零非对角线元素(行优先)\n" ); for( i = 0; i <= al_elem; i++ )
{
if( i == al_ln )
value = 0;
else
scanf( "%lf", &value );
matrix.data[i].value = value;
}
printf( "\n系数矩阵的信息相关如下(Quote为添加项,表示序列号):" );
printf( "\nV alue:" );
for( i = 0; i < al_elem + 1; i++ )
printf( "%3.0lf", matrix.data[i].value );
printf( "\nQuote:" );
for( i = 0; i < al_elem + 1; i++ )
printf( "%3d", i );
printf( "\n\n" );
for( i = 0; i <= al_elem; i++ )
{
if( i == 0 )
printf( "请参照上表输入第每行第一个非零非主对角元素的序号\n " );
if( i == al_ln )
matrix.data[i].index = al_elem + 1;
else
if( i == al_ln + 1 )
printf( "请依次输入非零非对角元素的列标\n " );
if( i != al_ln )
scanf( "%d", &matrix.data[i].index );
}
}
// 显示系数矩阵的相关信息
void Display_Matrix( struct Matrix matrix )
{
int i = 0;
printf( "系数矩阵采用MSR存储法的相关信息如下(Quote为添加项,表示序列号):" ); printf( "\nQuote: " );
for( i = 0; i <= al_elem; i++ )
printf( "%3d", i );
printf( "\nvalue: " );
for( i = 0; i <= al_elem; i++ )
printf( "%3.0lf", matrix.data[i].value );
printf( "\nIndex: " );
for( i = 0; i <= al_elem; i++ )
printf( "%3d", matrix.data[i].index );
printf( "\n" );
}
// 初始化向量,使之分量全为0
void Init_X( double *X, int n )
{
int i = 0;
for( i = 0; i < n; i++ )
X[i] = 0;
}
// 初始化右端向量
void Init_b( double *b, int n )
{
int i = 0;
Init_X( b, n );
printf( "\n请输入%d个数(方程组的右端项): ", n );
for( i = 0; i < n; i++ )
scanf( "%lf", &b[i] );
}
/
/ 初始化未知向量,得到初始值printf输出格式matlab
void Init_temp_X( double *temp, int n )
{
int i = 0;
Init_X( temp, n );
printf( "请输入X 向量的初始值: " );
for( i = 0; i < n; i++ )
scanf( "%lf", &temp[i] );
}
// 采用高斯-塞德尔GS迭代法求值
int Calculate_X_GS( struct Matrix A, double *X, double *b )
{
int i = 0, j = 0, sum = 0;
double temp = 0.0, temp_1 = 0.0, *temp_X = 0, judge = 0.0;
temp_X = (double *)new al_ln];
Init_temp_X( temp_X, A.total_ln );
for( i = 0; i < A.total_ln; i++ )
X[i] = temp_X[i];
sum = 0;
do
{
judge = 0.0;
for( i = 0; i < A.total_ln; i++ )
temp_X[i] = X[i];
for( i = 0; i < A.total_ln; i++ )
{
temp = 0.0;
temp_1 = 0.0;
for( j = A.data[i].index; A.data[j].index < i && j <= A.total_elem; j++ ) temp_1 += A.data[j].value*X[A.data[j].index];
for( ; j < A.data[i+1].index && j <= A.total_elem; j++ )
temp += A.data[j].value * temp_X[A.data[j].index];