数值分析计算实习题二.docx
- 文档编号:3243651
- 上传时间:2022-11-21
- 格式:DOCX
- 页数:18
- 大小:171.16KB
数值分析计算实习题二.docx
《数值分析计算实习题二.docx》由会员分享,可在线阅读,更多相关《数值分析计算实习题二.docx(18页珍藏版)》请在冰豆网上搜索。
数值分析计算实习题二
《数值分析》计算实习题二
算法设计方案
1.主要计算步骤:
计算函数f(x,y)在拟合所需的节点处的函数值。
将各拟合节点(xi,yj)分别带入非线性方程组
0.5cost+u+v+w–x=2.67
t+0.5sinu+v+w–y=1.07
0.5t+u+cosv+w–x=3.74
t+0.5u+v+sinw–y=0.79
解非线性方程组得解向量(tij,uij,vij,wij)。
对数表z(t,u)进行分片二次代数插值,求得对应(tij,uij)处的值,即为f(xi,yj)的值。
对上述拟合节点分别进行x,y最高次数为k(k=0,1,2,3…)次的多项式拟合。
每次拟合后验证误差大小,直到满足要求。
2.求解非线性方程组选择Newton迭代法,迭代过程中需要求解线性方程组,选择选主元的Doolittle分解法。
3.对z(t,u)进行插值选择分片二次插值。
4.拟合基函数φr(x)ψs(y)选择为φr(x)=xr,ψs(y)=ys。
拟合系数矩阵c通过连续两次解线性方程组求得。
一.源程序
#include"stdio.h"
#include"stdlib.h"
#include"math.h"
voidDoolittle(double*A,intn,int*M)
//功能说明:
对n阶矩阵A进行选主元的Doolittle分解
//参数说明:
A:
欲进行分解的方阵,同时也是返回参数,分解后的结果
//存储于A中
//n:
方阵A的维数
//M;(返回参数)n维向量,记录选主元过程中行交换的次序
{
inti,j,k,t;
double*s;
doubleMaxs,temp;
s=(double*)calloc(n,sizeof(double));
for(k=0;k { for(i=k;i { s[i]=A[i*n+k]; for(t=0;t } Maxs=abs(s[k]);M[k]=k; for(i=k+1;i { if(Maxs { Maxs=abs(s[i]); M[k]=i; } } if(M[k]! =k) { for(t=0;t { temp=A[k*n+t]; A[k*n+t]=A[M[k]*n+t]; A[M[k]*n+t]=temp; } temp=s[k]; s[k]=s[M[k]]; s[M[k]]=temp; } if(Maxs<(1e-14)) { s[k]=1e-14; printf("%.16e方阵奇异\n",Maxs); } A[k*n+k]=s[k]; for(j=k+1;(j { for(t=0;t A[j*n+k]=s[j]/A[k*n+k]; } } } voidSolve_LUEquation(double*A,intn,double*b,double*x) //功能说明: 解方程LUx=b,其中L、U共同存储在A中 //参数说明: A: 经Doolittle分解后的方阵 //n: 方阵A的维数 //b: 方程组的右端向量 //x: (返回参数)方程组的解向量 { inti,t; for(i=0;i { x[i]=b[i]; for(t=0;t } for(i=n-1;i>-1;i--) { for(t=i+1;t x[i]/=A[i*n+i]; } } voidTranspose(double*A,intm,intn,double*AT) //功能说明: 求m×n阶矩阵A的转置AT //参数说明: A: 已知m×n阶矩阵 //m: A的行数 //n: A的列数 //AT: (返回参数)A的转置矩阵(n×m) { inti,j; for(i=0;i for(j=0;j } voidSolve_LEquation(double*A,intn,double*B,double*x,intm) //功能说明: 解线性方程组Ax=B,该函数可对系数矩阵相同 //而右端向量不同的多个方程组同时求解。 //参数说明: A: 方程组系数矩阵 //n: A的维数 //B: m组右端向量构成的n×m矩阵 //x: (返回参数)方程组的解,n×m矩阵 //m: 不同右端向量的组数。 { int*M,i,j; M=(int*)calloc(n,sizeof(int)); double*BT,*xT,temp; BT=(double*)calloc(n*m,sizeof(double)); xT=(double*)calloc(n*m,sizeof(double)); //求B的转置BT,使得对应一个方程组的右端系数可以连续存储 //便于函数调用 Transpose(B,n,m,BT); Doolittle(A,n,M);//将A三角分解 for(i=0;i { for(j=0;j { temp=BT[i*n+j]; BT[i*n+j]=BT[i*n+M[j]]; BT[i*n+M[j]]=temp; } //对不同右端向量分别解方程组 Solve_LUEquation(A,n,BT+i*n,xT+i*n); } Transpose(xT,m,n,x); } voidMatrix_Mult(double*A,double*B,intm,ints,intn,double*C) //功能说明: 矩阵乘法AB=C,其中A为m×s阶,B为s×n阶。 //参数说明: A: m×s矩阵 //B: s×n矩阵 //m;A的行数 //s: A的列数,B的行数 //n: B的列数 //C: 返回值,m×n矩阵,返回AB的结果 { inti,j,k; for(i=0;i for(j=0;j { C[i*n+j]=0; for(k=0;k } } voidSet_NLtoL_JacobiA(double*A,double*x) //功能说明: 求题中非线性方程组对应自变量向量x的雅克比矩阵 //参数说明: A: 4×4矩阵 //x: 4维向量,非线性方程的自变量(t,u,v,w) { A[0]=-0.5*sin(x[0]);A[1]=1; A[2]=1;A[3]=1; A[4]=1;A[5]=0.5*cos(x[1]); A[6]=1;A[7]=1; A[8]=0.5;A[9]=1; A[10]=-sin(x[2]);A[11]=1; A[12]=1;A[13]=0.5; A[14]=1;A[15]=cos(x[3]); } voidSet_NLtoL_B(double*B,double*x,doublea,doubleb) //功能说明: 求非线性方程组Newton迭代法的右端式: -F(x) //参数说明: B: 4维向量,对应Newton迭代法的右端式-F(x) //x: 4维向量,非线性方程的自变量(t,u,v,w) //a: 非线性方程的参数,这里对应题中的x //b: 非线性方程的参数,这里对应题中的y { B[0]=-(0.5*cos(x[0])+x[1]+x[2]+x[3]-a-2.67); B[1]=-(x[0]+0.5*sin(x[1])+x[2]+x[3]-b-1.07); B[2]=-(0.5*x[0]+x[1]+cos(x[2])+x[3]-a-3.74); B[3]=-(x[0]+0.5*x[1]+x[2]+sin(x[3])-b-0.79); } doubleVector_FanShu(double*V,intn) //功能说明: 求n维向量V的无穷范数 //返回值: 所求范数 { inti; doublemax=0; for(i=0;i if(max returnmax; } voidSolve_NLEquation(doublea,doubleb,double*x) //功能说明: 求解非线性方程组 //参数说明: x: (返回值)4维向量,非线性方程的解向量(t,u,v,w) //a: 非线性方程的参数,这里对应题中的x //b: 非线性方程的参数,这里对应题中的y { doubleA[4][4],F[4],detx[4]; doubledet=1e-12;//求解精度要求 inti,k=0; intM=5000;//最大迭代次数 //设定迭代初始值 x[0]=0.5;x[1]=1;x[2]=1;x[3]=1; do { Set_NLtoL_B(F,x,a,b); Set_NLtoL_JacobiA(*A,x); Solve_LEquation(*A,4,F,detx,1); if(Vector_FanShu(detx,4)/Vector_FanShu(x,4) return; for(i=0;i<4;i++)x[i]+=detx[i]; k++; }while(k printf("Newton法在该初值不收敛\n"); } intNearPointIndex(double*v,doublea,intn) //功能说明: 查找n维向量V中与常数a最接近的元素的下标 //参数说明: V: n维向量 //n: 向量V的维数 //返回值: 与a最接近的元素的下标值 { inti,Index; doublemin; min=abs(v[0]-a); Index=0; for(i=1;i { if(min>abs(v[i]-a)) { min=abs(v[i]-a); Index=i; } } returnIndex; } doubleChaZhiZ(doublea,doubleb) //功能说明: 求数表z(t
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 计算 实习
![提示](https://static.bdocx.com/images/bang_tan.gif)