第1页 / 共19页
第2页 / 共19页
第3页 / 共19页
第4页 / 共19页
第5页 / 共19页
第6页 / 共19页
第7页 / 共19页
第8页 / 共19页
数值分析迭代法解非线性方程组c语言程序.doc
《数值分析 B》计算实习题目三
一、算法设计方案
,
(
i yxf
1. 确定数表
j
)
，其中 i=0,1,…,10；j=0,1,…,20 。
a) 将
(
,
i yx
j
)
代入非线性方程组解出
(
wvut
ij
ij
,
,
,
ij
ij
)
。
b) 利用所给的
z
,( uth
)
二维数表，用分片二次代数插值将上步得出的
,
(
ij ut
ij
)
值代入计算
(
,
ij uth
ij
)
的值。
c) 由
z
(
yxf
,
i
)
j
(
,
uth
ij
ij
)
得出数表
(
,
i yxf
j
)
2. 求
,(
yxf
)
在区域
D
,{(
yx
0|)
x
5.0,8.0
y
}5.1
上的一个近
似表达式
,(
yxp
)
k
,
sr
r
yxc
rs
0
s
，其中 r=0,1,…,k；s=0,1,…k 。
a) 以函数组
)(
x s
r
)(
y
r yx
s
为基函数，构成以{ rsc }为参数的曲面族
,(
yxp
)
k
,
sr
r
yxc
rs
0
s
进行曲面拟合。
b) 当曲面对数据的拟合精度
停止拟合。
10
20
0i
j
0
(
(
,
yxf
i
j
)
(
,
yxp
i
j
))
2
7
10
时
3. 求
(
,
yxf
*
i
*
j
),
(
,
yxp
*
i
*
j
)
的值,其中 i=1,2,…,8；j=1,2,…,5，以观察
,(
yxp
)
逼近
,(f
yx 的效果。
)
a) 用第 1 步中确定数表
(
,
i yxf
j
)
的方法将
x
*
i
.10
,
yi
*
j
05.05.0
j
代入求
(
xf
*
,
i y
*
j
)
。
b) 将
x
*
i
.10
,
yi
*
j
05.05.0
j
代入第 2 步求出的
,(
yxp
)
求
(
xp
*
,
i y
*
j
)
。
二、全部源程序
#include "stdio.h"
#include "math.h"
#include "stdlib.h"
#define E1 1e-12//---------------------------------------------------------------迭代求解精度
#define E2 1e-7//-----------------------------------------------------------------------拟合精度
int Gauss(double a[4][4], double b[4], double c[4]); //---------------------------函数声明
void Newton(double x, double y, double r[]);
double chazhi(double t, double u);
double pxy(double x, double y, int z);
int nihe();
int Doolittle(double a[11][11], int r);
void huidai1(double a[11][11], double b[11], double x[11], int r);
//------------------------------------------------------------------------------------------变量声明
double X[11], Y[21]; //-------------------------------------------------插值节点横, 纵坐标
double xyBiao[11][21];//------------------------------------------关于 z、x、y 的二维数表
double C[11][11];//---------------------------函数 z = f(x, y)的近似表达式的系数矩阵
int main()//---------------------------------------------------------------------------------主函数
{
double x, y, t, u;
double x1[8], y1[5], fxy[8][5], r[4],q[4];
int i, j;
for ( i = 0; i < 11; i ++)
X[i] = 0.08 * i;
for ( j = 0; j < 21; j ++)
Y[j] = 0.5 + 0.05 * j;
x = X[i];
for (j = 0; j < 21; j ++)
{
printf ("\n 数表:xi,yi,f(xi,yi)：\n");//-----------------建立关于 z、x、y 的二维数表
for (i = 0; i < 11; i ++)
{
y = Y[j];
q[0] = 0.5;
q[1] = 1.5;
q[2] = -0.3;
q[3] = 1.5;
Newton(x, y, q); //----------Newton 法解非线性方程组从 x，y 得到 t，u
t = q[0];
u = q[1];
xyBiao[i][j] = chazhi(t, u);
printf ("x =%1.2f ", X[i]);
printf (", y =%1.2f \n", Y[j]);
printf ("f(x,y)=%1.11e\n", xyBiao[i][j]);
}
}
int k = nihe();
//-------------计算 f(x[i], y[j]), p(x[i], y[j])的值，以观察 p(x, y)逼近 f(x, y)的效果
for (i = 0; i < 8; i ++)
x1[i] = 0.1 * (i + 1);
for (i = 0; i < 5; i ++)
y1[i] = 0.5 + 0.2 * (i + 1);
printf ("\n 数表 x*,y*,f(x*,y*), p(x*,y*)：\n");
for (i = 0; i < 8; i ++)
{
x = x1[i];
for (j = 0; j < 5; j ++)
{
}
}
return 0;
int i, j, k, r;
double m[4], w[4][5];
double t;
for (i = 0; i < 4; i ++)
{
for (j = 0; j < 4; j ++)
w[i][j] = a[i][j];
w[i][4] =b [i];
}
for (k = 0; k < 4; k ++)
{
y = y1[j];
r[0] = 0.5; //-------------------------------------------------------迭代初始向量
r[1] = 1.5;
r[2] = -0.3;
r[3] = 1.5;
Newton(x, y, r); //-----------Newton 法解非线性方程组从 x，y 得到 t，u
t = r[0];
u = r[1];
fxy[i][j] = chazhi(t, u);
printf ("x* =%1.1f ", x1[i]);
printf (", y*=%1.1f \n", y1[j]);
printf ("f(x*, y*) =%1.11e\n", fxy[i][j]);
printf ("p(x*, y*) =%1.11e\n", pxy(x, y, k));
}
//子函数一：----------------------------------------------------------列主元素 Gauss 消去法
int Gauss(double a[4][4], double b [4], double c[4])
{
r = k;
t = fabs(w[k][k]);
for (i = k + 1; i < 4; i ++)
{
if (fabs(w[i][k]) > t)
{
r = i;
t = fabs(w[i][k]);
}
}
if ( r != k)
{
for (j = k; j < 4 + 1; j ++)
{
t = w[k][j];
w[k][j] = w[r][j];
w[r][j] = t;
}
}
if (0 == w[k][k])
{
return 0;
}
for (i = k + 1; i < 4; i ++)
{
m[i] = w[i][k] / w[k][k];
for (j = k; j < 4 + 1; j ++)
{
w[i][j] = w[i][j] - m[i] * w[k][j];
}
}
}
for (i = 4 - 1; i >= 0; i --)
{
c [i] = w[i][4];
for (j = i + 1; j < 4; j ++)
c [i] -= w[i][j] * c[j];
c [i] /= w[i][i];
}
return 1;
}
//子函数二：------------------------------- Newton 法解非线性方程组从 x，y 得到 t，u
void Newton(double x, double y, double r[])
{
double c[4],s[4], a[4][4];
double e,max;
int i,j, d= 0,t=0,q=0;
do
{
d++;
for (i = 0; i < 4; i ++)//---------------------------------------------初始化系数矩阵
for (j = 0; j < 4; j ++)
a[i][j] = 1.0;
a[2][0] = a[3][1] = 0.5;
a[0][0] = -0.5 * sin(r[0]);
a[1][1] = 0.5 * cos(r[1]);
a[2][2] = -sin(r[2]);
a[3][3] = cos(r[3]);
c[0] = 2.67 - 0.5 * cos(r[0]) - r[1] - r[2] - r[3] + x;
c[1] = 1.07 - r[0] - 0.5 * sin(r[1]) - r[2] - r[3] + y;
c[2] = 3.74 - 0.5 * r[0] - r[1] - cos(r[2]) - r[3] + x;
c[3] = 0.79 - r[0] - 0.5 * r[1] - r[2] - sin(r[3]) + y;
if (!Gauss(a, c, s))
{
exit(-1);
}
max = fabs(s [0]);
for ( i = 1; i <4; i++)
{
if (fabs(s [i]) > max)
{
max = fabs(s [i]);
t= i;
}
}
max = fabs(r [0]);
for ( i = 1; i <4; i++)
{
if (fabs(r[i]) > max)
{
max = fabs(r [i]);
q= i;
}
}
e = fabs(s[t]) / fabs(r[q]);
for (i = 0; i < 4; i ++)
r[i] += s[i];
} while(e > E1);
double t1[3],u1[3];
int m,n, g, h;
m = int (10 * t) % 2 == 0 ? 5 * t : 5 * t + 1;
n = int (5 * t) % 2 == 0 ? 2.5 * t : 2.5 * t + 1;
m = 0 == m ? m + 1 : m;
m = 5 == m ? m - 1 : m;
n = 0 == n ? n + 1 : n;
n = 5 == n ? n - 1 : n;
for (g = m - 1; g < m + 2; g ++)
{
t1[g - m + 1] = 1.0;
for (h = m - 1; h < m + 2; h ++)
{
if (g != h)
{
}
}
t1[g - m + 1] *= t - t2[h];
t1[g - m + 1] /= t2[g] - t2[h];
}
//子函数三：---------------------------------------------------分片二次代数插值计算 h(t,u)
double chazhi(double t, double u)
{
=
{-0.50,
-0.34,
double
tuBiao[6][6]
2.06,
3.50,-0.42,-0.50,-0.26,0.30,1.18,2.38, -0.18, -0.50, -0.50, -0.18, 0.46,1.42, 0.22,
-0.34, -0.58, -0.50, -0.10,0.62, 0.78,-0.02, -0.50, -0.66, -0.50,-0.02, 1.50,0.46, -0.26,
-0.66, -0.74, -0.50 };//----------------------------------------------关于 z、t、u 的二维数表
double t2[6] = {0.0, 0.2, 0.4, 0.6, 0.8, 1.0};
double u2[6] = {0.0, 0.4, 0.8, 1.2, 1.6, 2.0};
0.14,
0.94,
}
for (g = n - 1; g < n + 2; g ++)
{
u1[g - n + 1] = 1.0;
for (h = n - 1; h < n + 2; h ++)
{
if (g != h)
{
u1[g - n + 1] *= u - u2[h];
u1[g - n + 1] /= u2[g] - u2[h];
}
}
}
double z = 0;
for (g = m - 1; g < m + 2; g ++)
for (h = n - 1; h < n + 2; h ++)
z += t1[g - m + 1] *u1[h - n + 1] * tuBiao[g][h];
return z;
}
//子函数四：------------------------------------------------计算 f(x, y)的近似表达式 p(x, y)
double pxy(double x, double y, int z)
{
double p = 0.0;
for(int r = 0; r <=z; r ++)
for (int s = 0; s <=z; s ++)
p += C[r][s] * pow(x, r) * pow(y, s);
return p;
}
//子函数五：----------------------------------------------------------------------------函数拟合
int nihe()
{
int i, j, k;
double wucha;
double A[11][21], B[11][11], V[11][11], D[21][11], G[21][11], Q[11][11];
printf ("\n 选择过程的 k 和 wucha 值：\n");
for (k = 1; k < 11; k ++)
{
{
}
}
for (i = 0; i < 11; i ++)
for (j = 0; j <= k; j ++)
{
B[i][j] = 1;
for (int t = 0; t < j; t ++)
B[i][j] *= X[i];
}
for (i = 0; i <= k; i ++)
{
for (j = 0; j <= k; j ++)
{
V[i][j] = 0;
for (int t = 0; t < 11; t ++)
V[i][j] += B[t][i] * B[t][j];
}
if (!Doolittle(V, k + 1))
{
exit(-1);
}
for (j = 0; j < 21; j ++)
{
double x[11],
b[11];
for (i = 0; i <= k; i ++)
{
b[i] = 0;
for (int t = 0; t < 11; t++)
b[i] += B[t][i] *xyBiao[t][j];
}
huidai1(V, b, x, k + 1);
for (i = 0; i <= k; i ++)
A[i][j] = x[i];
{
}
}
for (i = 0; i < 21; i ++)
{
for (j = 0; j <= k; j ++)
{
G[i][j] = 1;
for (int t = 0; t < j; t ++)
{
G[i][j] *= Y[i];
}
}
}
for (i = 0; i <= k; i ++)
{
for (j = 0; j <= k; j ++)
{
Q[i][j] = 0;
for (int t = 0; t < 21; t ++)
Q[i][j] += G[t][i] * G[t][j];
}
}
if (!Doolittle(Q, k + 1))
{
exit(-1);
}
for (j = 0; j < 21; j ++)
{
double x[11],
b[11];
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
