北京理工大学徐特立学院数值分析大作业上机实验
北航数值分析大作业 第一题 幂法与反幂法
数 值 分 析（B ） 大 作 业（一）姓名： 学号： 电话：1、算法设计：①求1λ、501λ和s λ的值：s λ：s λ表示矩阵的按模最小特征值，为求得s λ直接对待求矩阵A 应用反幂法即可。
1λ、501λ：若矩阵A 的特征值满足关系 1n λλ<<且1n λλ≠，要求1λ、及501λ时，可按如下方法求解： a ． 对矩阵A 用幂法，求得按模最大的特征值1m λ。
b ． 按平移量1m λ对矩阵A 进行原点平移得矩阵1m BA I λ=+，对矩阵B 用反幂法求得B 的按模最小特征值2m λ。
c ． 321m m m λλλ=-则：113min(,)m m λλλ=，13max(,)n m m λλλ=即为所求。
②求和A 的与数5011140k k λλμλ-=+最接近的特征值ik λ（k=0，1，…39）：求矩阵A 的特征值中与P 最接近的特征值的大小，采用原点平移的方法：先求矩阵 B=A-PI 对应的按模最小特征值k β，则k β+P 即为矩阵A 与P 最接近的特征值。
在本次计算实习中则是先求平移矩阵k B A I μ=-，对该矩阵应用反幂法求得s λ，则与k μ最接近的A 的特征值为：s P λ+重复以上过程39次即可求得ik λ（k=0，1，…39）的值。
③求A 的（谱范数）条件数2cond()A 和行列式det A ：在（1）中用反幂法求矩阵A 的按模最小特征值时，要用到Doolittle 分解方法，在Doolittle 分解完成后得到的两个矩阵分别为L 和U ，则A 的行列式可由U 阵求出，即：det(A)=det(U)。
求得det(A)不为0，因此A 为非奇异的实对称矩阵，则： max 2()scond A λλ=，max λ和s λ分别为模最大特征值与模最小特征值。
2、程序源代码：#include "Stdio.h"#include "Conio.h"#include "math.h"//****************************************************************************// // 在存储带状矩阵时，下面的几个量在程序中反复用到，为方便编程故把它们定义成宏.// // M ：转换后的矩阵的行数，M=R+S+1。
孙志忠北京理工大学偏微分方程数值解上机作业
偏微分方程数值解大作业目录第一题 (3)第二题 (7)第三题 (16)第四题 (20)第五题 (26)第六题（附加题1） (39)第七题（附加题2） (45)第八题（附加题3） (51)第一题习题13.（1）解曲线图图1 （2）误差曲线图图2（3）表格表1 部分点处精确解和取不同步长时所得的数值解表2 取不同步长时部分结点处数值解的误差的绝对值和数值解的最大误差（4）MATLAB源代码M=64;a=0;b=pi/2;h=(b-a)/M;x=[a+h:h:b-h];u=zeros(M-1,M-1);u(1,1)=(2/h^2)+(x(1)-1/2)^2;u(1,2)=-(1/h^2);u(M-1,M-1)=(2/h^2)+(x(M-1)-1/2)^2;u(M-1,M-2)=-(1/h^2);for i=2:M-2u(i,i-1)=-(1/h^2);u(i,i)=(2/h^2)+(x(i)-1/2)^2;u(i,i+1)=-(1/h^2);endf=zeros(M-1,1)f(1)=(x(1).*x(1)-x(1)+5/4).*sin(x(1));f(M-1)=(x(M-1).*x(M-1)-x(M-1)+5/4).*sin(x(M-1))+1/h^2; for j=2:M-2f(j)=(x(j).*x(j)-x(j)+5/4).*sin(x(j));endy=inv(u)*f; true=sin(x); plot(x,y'-true)第二题习题二（P67 第3题）（1）h=1/4, τ=1/4精确解数值解误差（2）h=1/8, τ=1/8精确解数值解误差（3）h=1/16, τ=1/16精确解数值解误差（4）h=1/32, τ=1/32精确解数值解误差（5）h=1/64, τ=1/64精确解精确解误差（6）表格（7）Matlab代码function[p,e,u,x,y,k]=fivepoint(h,m,n,kmax,ep)%五点差分法和G-S迭代法解椭圆型方程%kmax为最大迭代次数；%m,n分别为x,y方向的网格数；%ep为精度；%u为差分解，p为精确解，e为误差；%例如在命令行窗口输入[p,e,u,x,y,k]=fivepoint(1/64,64,64,10000,1e-10);%代表步长为1/64,精度为10^(-10),最大迭代次数为10000的五点差分格式syms temp;u=zeros(n+1,m+1);x=0+(0:m)*h;y=0+(0:n)*h;for(j=1:n+1)u(j,1)=sin(y(j))+cos(y(j));u(j,m+1)=exp(1)*(sin(y(j))+cos(y(j)));endfor(i=1:m+1)u(1,i)=exp(x(i));u(n+1,i)=exp(x(i))*(sin(1)+cos(1));endt=zeros(n-1,m-1);for(k=1:kmax)for(j=2:n)for(i=2:m)temp=(u(j,i+1)+u(j,i-1)+u(j+1,i)+u(j-1,i))/4; t(j,i)=(temp-u(j,i))*(temp-u(j,i));u(j,i)=temp;endendt(j,i)=sqrt(t(j,i));if(k>kmax)break;endif(max(max(t))<ep)break;endendfor(j=1:n+1)for(i=1:m+1)p(j,i)=(sin(y(j))+cos(y(j)))*exp(x(i));e(j,i)=abs(u(j,i)-p(j,i));endend代码使用说明：在命令行窗口输入[p,e,u,x,y,k]=fivepoint(1/64,64,64,10000,1e-10);代表步长为1/64,精度为10^(-10),最大迭代次数为10000的五点差分格式surf(x,y,p)可作出近似解曲线surf(x,y,u)可作出精确解曲线surf(x,y,e)可作出误差曲线注意精度不要选取的太低，否则随着步长减小误差反而增大。
北理工数据结构上机实验报告1
数据结构上级实验报告1——采用单向环表实现约瑟夫环一实验目的采用单向环表实现约瑟夫环二实验内容请按以下要求编程实现：①从键盘输入整数m，通过create函数生成一个具有m个结点的单向环表。
环表中的结点编号依次为1，2，……，m。
②从键盘输入整数s（1<=s<=m）和n，从环表的第s个结点开始计数为1，当计数到第n个结点时，输出该第n结点对应的编号，将该结点从环表中消除，从输出结点的下一个结点开始重新计数到n，这样，不断进行计数，不断进行输出，直到输出了这个环表的全部结点为止。
例如，m=10，s=3，n=4。
则输出序列为：6，10，4，9，5，2，1，3，8，7。
三程序设计(1)概要设计抽象数据类型：struct node{ int number;struct node * next;};宏：无主函数流程与各函数模块调用关系：创建头结点=〉输入链表的长度m=〉调用create函数建立链表=〉输入首次定位的节点s 以及间隔数n=〉调用deal函数处理链表=〉结束函数模块间无彼此调用关系(2)详细设计建立链表算法实现：void create(node * head,int m){node * p,*q;int i = 1;for(;m>0;m--){ p = (node *)malloc(sizeof(node));if(p == NULL) exit(0);p->number = m;p->next = head->next;head->next = p;if(i){q = p ; i--;}}q->next = head->next;}处理链表算法实现：void deal(node * head,int m,int s,int n){node * p,*q;int count=1,tag=0;p = head;while(s--){ p = p->next;}while(tag != m){if(count == n){printf("%d->",p->number);for(q = p->next; q->next != p; q=q->next);q->next = p->next;p = p->next;q = q->next;count = 1;tag++;}else {p = p->next;count++;}}printf("<-end\n");}主函数代码实现：int main(){ int m,s,n;node * head, * p;p = (node *)malloc(sizeof(node));p->next == NULL;head = p;printf("please input number m:\n");scanf("%d",&m);create(head,m);printf("please input number s(1<=s<=m):\n");scanf("%d",&s);printf("please input number n:\n");scanf("%d",&n);deal(head,m,s,n);system("pause");}四程序调试分析1在deal函数中进行对p节点的删除操作时出现问题，定位p的前一个节点q时，依旧采用了head指针从头遍历的方法，但是忽略了在函数操作中会删除掉头结点，这样head指针指向位置改变不能得到正确结果。
BUAA数值分析大作业三
北京航空航天大学2020届研究生《数值分析》实验作业第九题院系：xx学院学号：姓名：2020年11月Q9:方程组A.4一、 算法设计方案（一）总体思路1.题目要求∑∑===k i kj s r rsy x cy x p 00),(对f(x, y) 进行拟合，可选用乘积型最小二乘拟合。
),(i i y x 与),(i i y x f 的数表由方程组与表A-1得到。
2.),(**j i y x f 与1使用相同方法求得，),(**j i y x p 由计算得出的p(x,y)直接带入),(**j i y x 求得。
1. ),(i i y x 与),(i i y x f 的数表的获得对区域D ={ （x,y)|1≤x ≤1.24,1.0≤y ≤1.16}上的f (x , y )值可通过xi=1+0.008i ，yj=1+0.008j ，得到),(i i y x 共31×21组。
将每组带入A4方程组，即可获得五个二元函数组，通过简单牛顿迭代法求解这五个二元数组可获得z1~z5有关x,y 的表达式。
再将),(i i y x 分别带入z1~z5表达式即可获得f(x,y)值。
2.乘积型最小二乘曲面拟合2.1使用乘积型最小二乘拟合，根据k 值不用，有基函数矩阵如下：⎪⎪⎪⎭⎫ ⎝⎛=k i i k x x x x B 0000 ， ⎪⎪⎪⎭⎫ ⎝⎛=k j jk y y y y G 0000数表矩阵如下：⎪⎪⎪⎭⎫⎝⎛=),(),(),(),(0000j i i j y x f y x f y x f y x f U记C=[rs c ]，则系数rs c 的表达式矩阵为：11-)(-=G G UG B B B C T TT ）（通过求解如下线性方程，即可得到系数矩阵C 。
UG B G G C B B T T T =)()(2.2计算),(),,(****j i j i y x p y x f (i =1,2,…,31 ; j =1,2,…,21) 的值),(**j i y x f 的计算与),(j i y x f 相同。
数值分析上机实验
刘力辉 2010210804011.“画圆为方”问题也是古希腊人所提出几何三大难题中的另一个问题。
即求作一个正方形，使其面积等于已知圆的面积。
不妨设已知圆的半径为 R = 1，试用数值试验显示“画圆为方”问题计算过程中的误差。
（1）MATLAB 程序：y=pi^(1/2); % to generate 15-bit value of square root of pi b=1; d=1; for k=1:8 b=b*10;d=d/10; % b and d combined to control the digit of x x=d*fix(b*y); s(k)=x^3; l(k)=x; endformat long [l',s']（2）误差分析：2． 设，I x xd x n n=+⎰501由 x n = x n + 5 x n – 1 – 5 x n – 1 可得递推式 I n = – 5I n – 1 + 1/ n（1）从 I 0 尽可能精确的近似值出发，利用递推公式：I I nn n =-+-511 ( n = 1,2,…20)计算从 I 1 到 I 20 的近似值；（2）从 I 30 较粗糙的估计值出发，用递推公式：I I nn n -=-+11515 ( n= 30,29, …, 3, 2 )计算从I 1 到 I 20 的近似值；(3) 分析所得结果的可靠性以及出现这种现象的原因。
I 0 =dx x⎰+1051=ln (5+x )10|=ln 6-ln 5 所以I0≈0.18232155679395format longI0=log2(6)/log2(exp(1))-log2(5)/log2(exp(1)) % calculate the value of I0=ln6-ln5 for n=1:20I0=-5*I0+1/n; % recycling equation between I(n+1) and I(n) s(n)=I0; end s'则计算结果为：表1I1 0.0883922160302300 I11 0.0140713362538500 I2 0.0580389198488700 I12 0.0129766520640700 I3 0.0431387340890000 I13 0.0120398166027400 I4 0.0343063295550100 I14 0.0112294884148600 I5 0.0284683522249700 I15 0.0105192245923700 I6 0.0243249055418100 I16 0.0099038770381400 I7 0.0212326151481100 I17 0.0093041442210800 I8 0.0188369242594600 I18 0.0090348344501700 I9 0.0169264898137900 I19 0.0074574066965100 I10 0.0153675509310500I20 0.0127129665174600从计算的数据看出I 20=0.0127129665174600 > I 19=0.0074574066965100又I n 的积分范围为0～1，所以应该有I n >I n+1。
谷根代数值分析--上机实习报告
《数值分析》实验报告实验题目 非线性方程求解专业、班级 硕电力学号 姓名 课程编号 实验类型验证性 学时实验（上机）地点完成时间 任课教师谷根代 评分一、实验目的及要求1 掌握非线性方程数值解法的基本概念、方法；2 理解二分法、牛顿法、割线法、斯蒂芬森法求解非线性方程的原理、方法、步骤，会用计算机进行编程求解；3 掌握运用Matlab 进行数值计算的能力。
二、研究、解答以下问题问题：分别用二分法、牛顿迭代法、割线法、斯蒂芬森迭代法求方程26()(1)(1)0f x x x =+-=的根1x =，观察不同初试值条件下的收敛性，并给出你的结论【解】：2.1 算法分析 （1） 二分法二分法伪代码如下所示;输入：,a b ，定义函数()f x ，给定精度ε 输出：待解方程()f x 的根x 。
()/2x a b =+while (()0)&()f x n N ≠≤if ()()0f a f b <b x =()/21;else a x end x a b n n end==+=+由于本题中()f x 恒大于0，不满足二分法()()0f a f b <的条件，故本例无法通过二分法求出方程的根。
（2）牛顿迭代法牛顿迭代法迭代法迭代公式为：1'()()k k k k f x x x f x +=-a).对单根条件下，牛顿迭代法平方收敛，牛顿迭代法程序框图如下所示：b). 对重根条件下，此时迭代公式修改为1'(),()k k k k f x x x m m f x +=-⨯为根的重数 此时，牛顿迭代法至少平方收敛。
（3）割线法割线法的迭代公式为：111()(),1,2()()k k k k k k k f x x x x x k f x f x +--=--=-…]割线法是超线性收敛，其程序流程图为（4）斯蒂芬森迭代法斯蒂芬森迭代法的计算格式为21()()()2n n n n n n n n n n ny x z y y x x x z y x ϕϕ+==-=--+斯蒂芬森迭代法程序流程图与牛顿迭代法类似。
北京理工大学数值分析总复习
考试时带计算器
• 上机题请在11月30日晚9:30之前交，交 打印稿。
• 答疑时间：11月28，29, 30（即星期3, 4， 5）晚上7:30—9:30，上机作业也在答疑 时间交。
• 答疑地点：中教816。
2
第一章 误差
绝对（相对）误差 ( 限 ) 有效数字
(2) H H '((xx ii)) yyii,', (i0,1, ,n).
24
➢ H ( x ) y 0 h 0 ( x ) y 1 h 1 ( x ) y n h n ( x ) y 0 ' H 0 ( x ) y 1 ' H 1 ( x ) y n ' H n ( x )
(i 1 ,2 , ,n )
12
➢ 迭代法收敛的充分必要条件
x(k1) Mx(k) g,
x(0) 任意
收敛
(M)1.
➢ 迭代法收敛的充分条件
若i迭代法和Gauss-Seidel迭代法均收敛.
若A为对称正定阵, 则求解Ax=b的Gauss-Seidel迭代 法收敛.
复化梯形公式 复化Simpson公式 Romberg算法 Gauss型求积公式 代数精确度 截断误差
33
代数精确度
设有求积公式
b
n
f(x)dx
a
Akf(xk)
k0
若它对 f (x)=1, x, x2,…, xm 都能精确成立(即上式等
号成立), 但对 f (x)=xm+1 上式等号不成立, 则称该求
hi (x)
xi x0
(2n+1)次多项式
x i x n 1
h i ( x ) 1 2 l i ' ( x ) x ( x i ) l i 2 ( x ),
数值分析上机(C++版)
数值分析论文NUMERICAL ANALYSIS(THESIS)题目数值分析论文学生指导教师鸿雁学院专业班级2016年12月目录一二分法与牛顿迭代法1二拉格朗日插值法6三最小二乘法10四复合辛普森求积公式15五改进欧拉公式18六列主元消去法21七不动点迭代法26一二分法与牛顿迭代法一问题简介在有气体参加的恒容反应体系里，反应体系的总压会随着反应的进行而改变。
通过反应的平衡常数可以求出反应后气体的总量，而气体分子的分压与其所占比率是成正比的。
这类问题中关键是要求出反应后各分子的分压，各分子分压之和即为总压，其实质即为求方程的根。
例如：在298K下，化学反应2OF2=O2+2F2的平衡常数为0.410atm，如在298K下将OF2通入容器，当t=0时为1atm，问最后总压是多少？计算精度为10-3。
二数学模型假设是理想气体，可知2OF2=O2+2F2设氧的分压为p，平衡时有1-2p p p平衡时有整理得4p3-1640p2+164p-0410=0函数关系式为三算法选择与算法过程由计算得f(0.2)=－0.1156，f(0.3)=0.0424因此，有根区间为[0.2，0.3]用求单根的二分法计算时：可令a=0.2，b=0.3；最后计算得氧气分压：p=0.274609atm，总压为：P=3p+(1-2p)=1.274609atm用牛顿迭代法计算时：令则迭代公式为：令进行迭代，最后计算得氧气分压为：p=0.274901，则总压为：P=3p+(1-2p)=1.274901二分法的计算算法：#include <iostream>#include<math.h>using namespace std;double a=0.2,b=0.3,e=0.001,n;double f(float c){n=4*c*c*c-1.6400000*c*c+1.64*c-0.410;return n;}void main(){double c,m,p,a=0.2,b=0.3;for(p=a;b-a>=0.001;){c=(a+b)/2;m=f(c);if(m<0){a=c;}else if(m>0){b=c;}else{p=c;break;}p=(a+b)/2;}cout<<"\n a="<<a<<endl;cout<<"\n b="<<b<<endl;}牛顿迭代法的算法：#include<iostream>#include<string>#include<cmath>using namespace std;double f(double x){double m,n;m=4*x*x*x-1.640*x*x+1.64*x-0.410;n=12*x*x-3.280*x+1.64;return x-m/n;}void newton(double x,double d){double a=x;double b=f(a);int k; //记录循环的次数for(k=1;fabs(a-b)>d;k++){a=b;b=f(a);if(k>100){cout<<"迭代失败，该函数可能不收敛!"<<endl;}}cout<<"\n b="<<b<<endl;cout<<"\n k="<<k<<endl;return;}int main(){cout<<"请输入初始值x0和要求得结果的精度：";double x,d;cin>>x>>d;newton(x,d);return 0;}四数值实验过程二分法通过VC6.0程序运行的结果如下图所示：牛顿迭代法通过VC6.0的计算结果如下：五相关数值分析和实际应用分析由二分法程序输出结果可得:a=0.274219,b=0.275。
数值分析上机作业
数值分析上机作业2实验一：（1） ①用不动点迭代法求()013=--=x x x f 的根，至少设计两种迭代格式，一个收敛一个发散，1210-=ε。
（2） ②对迭代格式使用Aitken 加速，观察敛散性变化。
1取递推公式31)1(+=x x ，可以得到收敛时的迭代结果为：x=(2)^(1/3); t=1; while(1)if(abs(x-(x+1)^(1/3))<10^-12) break; endx=(x+1)^(1/3);t=t+1;end t xtemp=x^3-x-1 %带回来验证下 t = 16 x =1.324717957243755 temp =-4.225952920933196e-12 加速后代码如下 x=1;x=(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)) t=1; while(1)if(abs(x-(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)))<10^-10) break; endx=(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)); t=t+1; if (t>100000)break; %防止运算次数过多 end endfprintf('需要%d 次',t);输出需要670次>>此处取10^10-=ε，若到10^-12次方则可能需要运行更多次取13-=x x 则迭代发散。
使用aitken 加速计算结果如下 x=2; t=1; while(1)if( abs(x-((x*((x^3-1)^3-1))-(x^3-1)^2)/(x-2*(x^3-1)+(x^3-1)^3-1))<10^-10) break; end t=t+1;x=((x*((x^3-1)^3-1))-(x^3-1)^2)/(x-2*(x^3-1)+(x^3-1)^3-1); if (t>100000)break; %防止运算次数过多 end end t xt = 108x = 1.324717956244172由此可见经过aitken 加速以后，原来发散的迭代格式收敛了。
北理工数值分析大作业
数值分析上机作业第 1 章1.1计算积分，n=9。
（要求计算结果具有6位有效数字）程序：n=1:19;I=zeros(1,19);I(19)=1/2*((exp(-1)/20)+(1/20));I(18)=1/2*((exp(-1)/19)+(1/19));for i=2:10I(19-i)=1/(20-i)*(1-I(20-i));endformat longdisp(I(1:19))结果截图及分析：在MATLAB中运行以上代码，得到结果如下图所示：当计算到数列的第10项时，所得的结果即为n=9时的准确积分值。
取6位有效数字可得.1.2分别将区间[-10.10]分为100,200,400等份，利用mesh或surf命令画出二元函数z=的三维图形。
程序：>> x = -10:0.1:10;y = -10:0.1:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长0.1')>> x = -10:0.2:10;y = -10:0.2:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长 0.2')>>x = -10:0.05:10;y = -10:0.05:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长0.05')结果截图及分析：由图可知，步长越小时，绘得的图形越精确。
数值分析上机实验报告
数值分析上机实验报告实验报告插值法与数值积分实验（数值计算方法，3学时）一实验目的1．掌握不等距节点下的牛顿插值公式以及拉格朗日插值公式。
2．掌握复化的梯形公式、辛扑生公式、牛顿-柯特斯公式计算积分。
3. 会用龙贝格公式和高斯公式计算积分。
二实验内容用拉格朗日插值公式计算01.54.1==y x 以及所对应的近似值。
用牛顿插值公式求)102(y 的近似值。
三实验步骤（算法）与结果1拉格朗日插值法：（C 语言版）#include "Stdio.h" #include "Conio.h"int main(void) {float X[20],Y[20],x; int n;void input(float *,float *,float *,int *); float F(float *,float *,float,int); input(X,Y,&x,&n);printf("F(%f)=%f",x,F(X,Y,x,n));getch(); return 0; }void input(float *X,float *Y,float *x,int *n) {int i;printf("Please input the number of the data:");scanf("%d",n);printf("\nPlease input the locate of each num:\n");for(i=0;i<*n;i++){scanf("%f,%f",X+i,Y+i);}printf("\nPlease input the chazhi:"); scanf("%f",x);}float F(float *X,float *Y,float x,int n){int i,j;float Lx,Fx=0;for(i=0;i<n;i++)< p="">{Lx=1;for(j=0;j<n;j++)< p="">{if(j!=i) Lx=Lx*((x-*(X+j))/(*(X+i)-*(X+j))); } Fx=Fx+Lx*(*(Y+i));}return Fx;}得出结果如图：所以Y(1.4)=3.7295252#include#define N 10double X[N], Y[N], A[N][N];int n;double Newton(double x);double f(double x);void main() {printf("请输入已知x与对应y=f(x)的个数: n = "); scanf("%d", &n);getchar();if(n>N||n<=0) {printf("由于该维数过于犀利, 导致程序退出!"); return;}printf("\n请输入X[%d]: ", n);for (int i=0; i<="" p="">scanf("%lf", &X[i]);getchar();printf("\n请输入Y[%d]: ", n);for (i=0; i<="" p="">scanf("%lf", &Y[i]);getchar();double x;printf("\n请输入所求结点坐标x = ");scanf("%lf", &x);getchar();printf("\nf(%.4lf)≈%lf\n\n", x, Newton(x));}double Newton(double x) {int i, j;// 求均差for (i =0; i<="" p="">A[i][0] = Y[i];for (i=1; i<="" p="">for (j =1; j<=i; j++)A[i][j] = (A[i][j-1] - A[i-1][j-1]) / (X[i] - X[i-j]); // 求结点double result = A[0][0];for (i=1; i<="">double tmp = 1.0;for (int j=0; j<="" p="">tmp *= (x - X[j]);result += tmp * A[i][i];}return result;}四实验收获与教师评语</n;j++)<></n;i++)<>。
数值计算方法丁丽娟课后习题答案
数值计算方法丁丽娟课后习题答案【篇一：北京理工大学数值计算方法大作业数值实验1】）书p14/4分别将区间[?10,10]分为100,200,400等份，利用mesh或surf命令画出二元函数的三维图形。
z=|??|+ ??+?? +??++??【matlab求解】[x,y]=meshgrid(-10:0.1:10);a=exp(-abs(x));b=cos(x+y);c=1./(x.^2+y.^2+1);z=a+b+c;mesh(x,y,z);[x,y]=meshgrid(-10:0.05:10);a=exp(-abs(x));b=cos(x+y);c=1./(x.^2+y.^2+1);z=a+b+c;mesh(x,y,z);[x,y]=meshgrid(-10:0.025:10); a=exp(-abs(x));b=cos(x+y);c=1./(x.^2+y.^2+1);z=a+b+c;mesh(x,y,z);（二）书p7/1.3.2数值计算的稳定性（i）取= ??c语言程序—不稳定解 +=ln1.2，按公式=?? （n=1,2，…） #includestdio.h#includeconio.h#includemath.hvoid main(){float m=log(6.0)-log(5.0),n;int i;i=1;printf(y[0]=%-20f,m); while(i20){n=1/i-5*m;printf(y[%d]=%-20f,i,n);m=n;i++;if (i%3==0) printf(\n); }getch();}（ii） c语言程序—稳定解≈??[ ??+?? +?? ??+??按公式 =??(??)#includestdio.h#includeconio.h#includemath.hvoid main(){float m=(1/105.0+1/126.0)/2,n; k=n,n-1,n-2，…）（【篇二：北京理工大学数值计算方法大作业数值实验4】 p260/1考纽螺线的形状像钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为= ?????????????????? ?? ??????????= ?????????????? ??曲线关于原点对称，取a=1，参数s的变化范围[-5,5]，容许误差限分别是，,和。
北理工数值分析大作业
北理工数值分析大作业近年来，随着计算机技术的飞速发展，数值分析成为了计算机科学与工程中的重要领域。
数值分析着重于通过数值方法和计算技术解决实际问题，并通过数值计算的精度、稳定性和效率来评价算法的优劣。
本文将介绍北理工数值分析大作业的主题和方案。
本次数值分析大作业的主题是解决常微分方程（Ordinary Differential Equation，ODE）的初值问题。
常微分方程是描述自然界中许多物理现象和工程问题的数学模型，解决常微分方程的初值问题是数值分析中的经典问题之一我们将采用龙格-库塔法（Runge-Kutta method）来解决常微分方程的初值问题。
龙格-库塔法是一种常用的数值求解ODE的方法，其核心思想是通过将微分方程近似为差分方程，用离散化的方法逼近真实解。
具体方法是以一种递归的方式计算不同阶的近似解，从而提高解的精度。
首先，我们将介绍龙格-库塔法的基本原理和算法。
龙格-库塔法的核心是通过递归计算中间值来逼近真实解，并根据这些中间值计算最终的近似解。
我们将详细介绍四阶龙格-库塔法的具体算法和推导过程，并给出其误差分析和收敛性证明。
接下来，我们将选取一个具体的常微分方程作为案例来验证龙格-库塔法的有效性和稳定性。
我们选择经典的二阶线性常微分方程作为案例，该方程具有已知的解析解。
我们将通过比较数值解和解析解的差异来评估龙格-库塔法的精度和可靠性，并对算法的参数设置进行优化。
最后，我们将对龙格-库塔法进行性能测试和分析。
我们将使用不同的问题规模和不同的算法参数来测试算法的运行时间和内存占用情况，并通过绘制性能曲线来评估算法的效率和可扩展性。
总结起来，本文介绍了北理工数值分析大作业的主题和方案，主要包括龙格-库塔法的基本原理和算法、常微分方程的解析解验证、算法参数优化以及性能测试。
通过本次大作业的完成，我们将对数值分析的基本理论和方法有更深入的理解，并能够熟练应用数值方法解决实际问题。
数值分析上机实验报告
一、实验目的通过本次上机实验，掌握数值分析中常用的算法，如二分法、牛顿法、不动点迭代法、弦截法等，并能够运用这些算法解决实际问题。
同时，提高编程能力，加深对数值分析理论知识的理解。
二、实验环境1. 操作系统：Windows 102. 编程语言：MATLAB3. 实验工具：MATLAB数值分析工具箱三、实验内容1. 二分法求方程根二分法是一种常用的求方程根的方法，适用于连续函数。
其基本思想是：从区间[a, b]中选取中点c，判断f(c)的符号，若f(c)与f(a)同号，则新的区间为[a, c]，否则为[c, b]。
重复此过程，直至满足精度要求。
2. 牛顿法求方程根牛顿法是一种迭代法，适用于可导函数。
其基本思想是：利用函数在某点的导数值，求出函数在该点的切线方程，切线与x轴的交点即为方程的近似根。
3. 不动点迭代法求方程根不动点迭代法是一种迭代法，适用于具有不动点的函数。
其基本思想是：从初始值x0开始，不断迭代函数g(x)的值，直至满足精度要求。
4. 弦截法求方程根弦截法是一种线性近似方法，适用于可导函数。
其基本思想是：利用两点间的直线近似代替曲线，求出直线与x轴的交点作为方程的近似根。
四、实验步骤1. 二分法求方程根（1）编写二分法函数：function [root, error] = bisection(a, b, tol)（2）输入初始区间[a, b]和精度要求tol（3）调用函数计算根：[root, error] = bisection(a, b, tol)2. 牛顿法求方程根（1）编写牛顿法函数：function [root, error] = newton(f, df, x0, tol)（2）输入函数f、导数df、初始值x0和精度要求tol（3）调用函数计算根：[root, error] = newton(f, df, x0, tol)3. 不动点迭代法求方程根（1）编写不动点迭代法函数：function [root, error] = fixed_point(g, x0, tol)（2）输入函数g、初始值x0和精度要求tol（3）调用函数计算根：[root, error] = fixed_point(g, x0, tol)4. 弦截法求方程根（1）编写弦截法函数：function [root, error] = secant(f, x0, x1, tol)（2）输入函数f、初始值x0和x1，以及精度要求tol（3）调用函数计算根：[root, error] = secant(f, x0, x1, tol)五、实验结果与分析1. 二分法求方程根以方程f(x) = x^2 - 2 = 0为例，输入初始区间[a, b]为[1, 3]，精度要求tol 为1e-6。
北理工微机原理上机实验考试题目类型及参考程序
北理工微机原理上机实验考试题目类型及参考程序1到100奇数的累加和显示：DATA SEGMENTSUM DB 4 DUP(?),'$'DATA ENDSCODE SEGMENTASSUME CS:CODE,DS:DATA START:MOV AX,DATA MOV DS,AXMOV AX,0MOV BX,1MOV CX,50AGAIN:ADD AX,BXADD BX,2LOOP AGAINMOV CX,4LEA SI,SUMMOV [SI],AXNEXT: MOV AX,[SI]ROL AX,4MOV [SI],AXAND AX,000FHCMP AL,9JA NEXT1ADD AL,30HMOV DL,ALMOV AH,2INT 21HJMP GOONNEXT1:ADD AL,37HMOV DL,ALMOV AH,2INT 21HGOON:LOOP NEXTMOV DL,48HMOV AH,2INT 21HMOV AH,4CHINT 21HCODE ENDSEND START八个数判断奇偶个数：DATA SEGMENTJI DB 0OU DB 0BUFFER DB 01H,02H,03H,04H,05H,06H,07H,08H DATA ENDS CODE SEGMENTASSUME CS:CODE,DS:DATASTART:MOV AX,DATAMOV DS,AXMOV CX,8LEA SI,BUFFERAGAIN:MOV AL,[SI]MOV BL,2MOV AH,00HDIV BLAND AH,AHJNZ JISHUINC OUJMP GONEJISHU:INC JIGONE: INC SILOOP AGAINMOV DL,JIADD DL,30HMOV AH,2INT 21HMOV DL,' 'MOV AH,2INT 21HMOV DL,OUADD DL,30HMOV AH,2INT 21HMOV AH,4CHINT 21HCODE ENDSEND START十个有符号数从大到小排序并显示DATA SEGMENTTABLE1 DB 11H,33H,22H,44H,55H,66H,77H,88H,99H,0AH DATA ENDSCODE SEGMENTASSUME CS:CODE,DS:DATASTART:MOV AX,DATAMOV DS,AXMOV CH,9NEXT:XOR AX,AX LEA SI,TABLE1 MOV CL,9 CLDNEXT1: LODSBCMP AL,[SI] JG NEXT2 XCHG AL,[SI] MOV [SI-1],AL NEXT2:DEC CLJNZ NEXT1 DEC CHJNZ NEXT LEA SI,TABLE1 MOV CX,10 NEXT3:XOR AX,AX CLD LODSBMOV BX,16 XOR DX,DX DIV BX XCHG DL,AL PUSH AX ADD DL,30H CMP DL,'9' JBE NEXT4 ADD DL,07HMOV AH,2INT 21HPOP AXMOV DL,ALADD DL,30HCMP DL,'9'JBE NEXT5ADD DL,07HNEXT5:MOV AH,2INT 21HMOV DL,'H'MOV AH,2INT 21HMOV DL,' 'MOV AH,2INT 21HLOOP NEXT3NEXT6:MOV AH,4CHINT 21HCODE ENDSEND START五个有符号数从小到大排列并显示DATA SEGMENTTABLE1 DB 11H,33H,22H,44H,55H,66H,77H,88H,99H,0AH DATA ENDSCODE SEGMENTASSUME CS:CODE,DS:DATAMOV AX,DATAMOV DS,AXMOV CH,4红色字体是和上一个程序的不同点NEXT: XOR AX,AXLEA SI,TABLE1MOV CL,4CLDNEXT1:LODSBCMP AL,[SI]JG NEXT2XCHG AL,[SI]MOV [SI-1],ALNEXT2:DEC CLJNZ NEXT1DEC CHJNZ NEXTLEA SI,TABLE1MOV CX,5 NEXT3:XOR AX,AXCLDLODSBMOV BX,16XOR DX,DXDIV BXXCHG DL,ALPUSH AXADD DL,30HCMP DL,'9'JBE NEXT4ADD DL,07H NEXT4:MOV AH,2INT 21HPOP AXMOV DL,ALADD DL,30HCMP DL,'9'JBE NEXT5ADD DL,07H NEXT5:MOV AH,2INT 21HMOV DL,'H'MOV AH,2INT 21HMOV DL,' 'MOV AH,2INT 21HLOOP NEXT3 NEXT6:MOV AH,4CHINT 21H CODE ENDSEND START输出最大值以ASC11码输出DATA SEGMENTMEM DB 12H,56H,89H,0ABH,0DFH,29H,0,87H RESULT DB ? DATA ENDSCODE SEGMENTASSUME CS:CODE,DS:DATASTART:MOV AX,DATAMOV DS,AXLEA SI,MEMMOV CX,8MOV AL,0 NEXT: CMP AL,[SI] JA L1MOV AL,[SI]L1: INC SILOOP NEXTMOV CX,2LEA SI,RESULT MOV [SI],AL DISP: MOV AL,[SI] ROL AL,4MOV [SI],ALAND AL,0FHCMP AL,9JA NEXT1ADD AL,30H MOV DL,ALMOV AH,2INT 21HJMP NEXT2 NEXT1:ADD AL,37H MOV DL,ALMOV AH,2INT 21HNEXT2:LOOP DISP MOV AH,4CHINT 21HEND START字符串中0,1的个数DATA SEGMENTNUM DW 3FFFHONE DB 0ZERO DB 0DATA ENDSCODE SEGMENTASSUME CS:CODE,DS:DATA START: MOV AX,DATAMOV DS,AXLEA SI,NUMLEA DI,ZEROMOV AX,[SI]LEA SI,ONEMOV CH,16NEXT1:ROL AX,1JC NEXT2INC BYTE PTR [DI]JMP NEXT3NEXT2:INC BYTE PTR [SI]NEXT3:DEC CHJNZ NEXT1NEXT4:MOV AH,4CHINT 21HEND START字符串中0,1的个数带输出DATA SEGMENT NUM DW 3FFFHONE DB 0ZERO DB 0DATA ENDSCODE SEGMENTASSUME CS:CODE,DS:DATA START: MOV AX,DATAMOV DS,AXLEA SI,NUMLEA DI,ZEROMOV AX,[SI]LEA SI,ONEMOV CH,16NEXT1:ROL AX,1JC NEXT2INC BYTE PTR [DI]JMP NEXT3NEXT2:INC BYTE PTR [SI]NEXT3:DEC CHJNZ NEXT1NEXT4:MOV BX,10XOR DX,DXMOV DL,[DI]ADD DL,30HMOV AH,2INT 21HMOV DL,' 'MOV AH,2INT 21HXOR DX,DXXOR AX,AXMOV AX,[SI]XOR AH,AHDIV BXXCHG DL,ALADD DL,30H PUSH AXMOV AH,2INT 21HPOP AXXOR DX,DXMOV DL,ALADD DL,30HMOV AH,2INT 21H NEXT5: MOV AH,4CHINT 21H CODE ENDS END START。
数值分析上机实验报告
end
%%%%计算消元,得三角阵
for r=k+1:n;
m=A(r,k)/A(k,k);
for q=k:n+1;
A(r,q)=A(r,q)-A(k,q)*m;
end
end
end
%%%%求解x
x(n)=A(n,n+1)/A(n,n);
for k=n-1:-1:1;
三次样条插值：
给定区间[a, b]一个分划
⊿：a=x0<x1<…<xN=b
若函数S(x)满足下列条件：
1)S(x)在每个区间[xi, xj]上是不高于3次的多项式。
2)S(x)及其2阶导数在[a, b]上连续。则称S(x)使关于分划⊿的三次样条函数。
程序设计：
本实验采用Matlab的M文件编写。其中待插值的方程写成function的方式，如下
for k=1:n+1;
if k~=q;
l=l.*(x-s(k))./(s(q)-s(k));
else
l=l;
end
end
f=f+Rf(s(q))*l;%求插值函数
end
plot(x,f,'r')%作出插值函数曲线
grid on
hold on
分段线性插值源程序
clear
n=input('将区间分为的等份数输入：\n');
并将第r行和第k行的元素进行交换，以使得当前的 的数值比0要大的多。这种列主元的消去法的主要步骤如下：
1.消元过程
对k=1,2,…,n-1,进行如下步骤。
1)选主元，记
若 很小，这说明方程的系数矩阵严重病态，给出警告，提示结果可能不对。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
北京理工大学徐特立学院数值分析课后上机实验选做教材：数值计算方法（2011第一版）.丁丽娟,程杞元.高等教育出版社^以下代码作者原创^超链接：1.2（题目，原理，截图，代码）2.2（题目，原理，截图，代码）3.1（题目，原理，截图，c代码，m代码）5.1（题目，原理，截图，代码）5.3（题目，原理，截图，代码）第一章：数值计算中的误差2、题目简介：利用pi/4=1-1/3+1/5-1/7。
级数计算pi的近似值。
输入：误差值输出：求和项数，并输出pi值工具：C语言运行环境：VC-6.0计算公式及原理：利用pi/4=1-1/3+1/5-1/7。
级数计算pi的近似值，由数学原理可知误差会小于首次舍弃的项,可以编写循环实现。
程序运行结果截图：程序代码：（c语言）#include<stdio.h>void main(){printf("第一章第2题求pi,欢迎使用，请按提示操作。
\n");int i=1,n=0,k=1;double e,pi,er;printf("请输入误差(例如1e-4):");scanf("%lf",&e);printf("请稍候。
\n");er=e;pi=0;while(er>=e){pi+=k*1.0/i;k=-k;er=1.0/i;i+=2;n++;}pi*=4;printf("%d项求和后可以达到%.10lf精度，这时pi=%.10lf\n",n,e,pi);getchar();getchar();}第二章：解线性方程组的直接方法2、题目简介：用MATLAB软件编程实现追赶法求解三对角方程组的算法，并考虑梯形电阻电路问题，电路如下：其中电路中的各个电流{1i ，2i ，…，8i }须满足下列线性方程组：R V i i =- 22 21 0 252321=-+-i i i 0 252 432=-+-i i i 0 252 543=-+-i i i 0 252 654=-+-i i i 0 252 765=-+-i i i 0 252 876=-+-i i i 052 87=+-i i设V 220=V ，Ω=27R ，运用求各段电路的电流量。
输入：三列对角线元素向量，右侧常数元素向量 输出：三对角方程的解 工具：m 语言运行环境：MATLAB R2012.b 计算公式与原理：1481.827220≈=R V 上述方程组可用矩阵表示为：⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--------------00000001481.8522520000002520000002520000002520000002520000002520000002287654321i i i i i i i i 根据三对角方程追赶法原理以及公式可解。
程序运行结果截图：程序代码：（matlab ）%说明：追赶法是适用于三对角矩阵的线性方程组求解的方法，并不适用于其他类型矩阵。
%定义三对角矩阵A 的各组成单元。
方程为Ax=d function x=zhuiganfaa=input('输入左下角对角元素a （n ）(例如[1 2 3]):'); b=input('输入对角元素b （n ）(例如[1 2 3 4])：');c=input('输入左上角对角元素c （n ）(例如[1 2 3])：');d=input('输入列矩阵d （n ）例如[1 2 3 4]：'); n=length(b); u0=0;y0=0;%“追”的过程 L(1)=b(1);y(1)=(d(1))/L(1); u(1)=c(1)/L(1); for i=2:(n-1)L(i)=b(i)-a(i-1)*u(i-1);y(i)=(d(i)-y(i-1)*a(i-1))/L(i); u(i)=c(i)/L(i); endL(n)=b(n)-a(n-1)*u(n-1);y(n)=(d(n)-y(n-1)*a(n-1))/L(n); %“赶”的过程 x(n)=y(n); for i=(n-1):-1:1x(i)=y(i)-u(i)*x(i+1); end第三章：解线性方程组的迭代法1、题目简介：试分别用（1）Jacobi 迭代法；（2）Gauss-Seidel 迭代法；（3）共轭梯度法解线性方程组⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡--=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡--------1217142712151534112323537123219143211054321x x x x x迭代初始向量取)0(x =(0，0，0，0，0T )。
输入：根据要求输入精度或者步数，系数矩阵，向量b ，初值向量等 输出：解的精度，方程的解，要达到精度所需步数 工具：C 语言/m 语言 运行环境：VC-6.0/MATLAB计算公式及原理：根据相关迭代原理和算法，设计程序求解，程序可以按步数测试，也可以按精度测试。
运行结果截图：(C语言)运行结果截图（MATLAB）程序代码:（c语言）Jacobi迭代法解线性方程组:#include<stdio.h>#include<math.h>void fun(double a[][100],double *b,double *x,double *y,double *max, int N) {int i,j;double sum;for(i=1;i<=N;i++)x[i]=y[i];for(i=1;i<=N;i++){sum=0;for(j=1;j<=N;j++){if(j!=i)sum+=(a[i][j])*(x[j]);}y[i]=(b[i]-sum)/a[i][i];}*max=fabs(y[1]-x[1]);for(i=1;i<=N;i++)if(fabs(y[i]-x[i])>*max)}void main(){printf("欢迎来到Jacobi迭代法解线性方程组，请先确保该法收敛，按提示操作*^_^*\n");int N,i,j,k=0,mode,num;double a[100][100],b[100],x[100],y[100],e,sum,max;printf("请输入结束条件，迭代次数限制请按1，误差限制请按2:");scanf("%d",&mode);if(mode==2){printf("请输入误差限e:");scanf("%lf",&e);}else{printf("请输入迭代次数：");scanf("%d",&num);}printf("请输入方程阶数：");scanf("%d",&N);printf("请输入系数矩阵：\n");for(i=1;i<=N;i++)for(j=1;j<=N;j++)scanf("%lf",&a[i][j]);printf("请输入向量b:\n");for(i=1;i<=N;i++)scanf("%lf",&b[i]);printf("请输入初始解x:");for(j=1;j<=N;j++)scanf("%lf",&x[j]);for(i=1;i<=N;i++){sum=0;for(j=1;j<=N;j++){if(j!=i)sum+=a[i][j]*x[j];}y[i]=(b[i]-sum)/a[i][i];}max=fabs(y[1]-x[1]);for(i=1;i<=N;i++)if(fabs(y[i]-x[i])>max)printf("%d ",k);for(i=1;i<=N;i++)printf("x[%d]=%-10.7lf",i,x[i]);printf("\n");k++;printf("%d ",k);for(i=1;i<=N;i++)printf("x[%d]=%-10.7lf",i,y[i]);printf("\n");if(mode==1){while(k<num){fun(a,b,x,y,&max,N);k++;if(k<10)printf("%d ",k);elseprintf("%d ",k);for(i=1;i<=N;i++)printf("x[%d]=%-10.7lf",i,y[i]);printf("\n");}}else{while(max>e){fun(a,b,x,y,&max,N);k++;if(k<10)printf("%d ",k);elseprintf("%d ",k);for(i=1;i<=N;i++)printf("x[%d]=%-10.7lf",i,y[i]);printf("\n");}}getchar();getchar();}Gauss-Sdidel迭代法解线性方程组:#include<stdio.h>#include<math.h>void fun(double a[][100],double *b,double *x,double *y,double *max, int N){int i,j;double sum,ee[100];for(i=1;i<=N;i++)x[i]=y[i];for(i=1;i<=N;i++){sum=0;for(j=1;j<=N;j++){if(j!=i)sum+=(a[i][j])*(x[j]);}y[i]=(b[i]-sum)/a[i][i];ee[i]=y[i]-x[i];x[i]=y[i];}*max=fabs(ee[1]);for(i=1;i<=N;i++)if(fabs(ee[i])>*max)*max=fabs(ee[i]);}void main(){printf("欢迎来到Gauss-Sdidel迭代法解线性方程组，请先确保该法收敛，按提示操作*^_^*\n");int N,i,j,k=0,mode,num;double a[100][100],b[100],x0[100],x[100],y[100],e,ee[100],sum,max;printf("请输入结束条件，迭代次数限制请按1，误差限制请按2:");scanf("%d",&mode);if(mode==2){printf("请输入误差限e:");scanf("%lf",&e);}else{printf("请输入迭代次数：");scanf("%d",&num);}printf("请输入方程阶数：");scanf("%d",&N);printf("请输入系数矩阵：\n");for(i=1;i<=N;i++)for(j=1;j<=N;j++)scanf("%lf",&a[i][j]);printf("请输入向量b:\n");for(i=1;i<=N;i++)scanf("%lf",&b[i]);printf("请输入初始解x:");for(j=1;j<=N;j++){scanf("%lf",&x[j]);x0[j]=x[j];}for(i=1;i<=N;i++){sum=0;for(j=1;j<=N;j++){if(j!=i)sum+=a[i][j]*x[j];}y[i]=(b[i]-sum)/a[i][i];ee[i]=y[i]-x[i];x[i]=y[i];}max=fabs(ee[1]);for(i=1;i<=N;i++)if(fabs(ee[i])>max)max=fabs(ee[i]);printf("%d ",k);for(i=1;i<=N;i++)printf("x[%d]=%-10.7lf",i,x0[i]); printf("\n");k++;printf("%d ",k);for(i=1;i<=N;i++)printf("x[%d]=%-10.7lf",i,x[i]); printf("\n");if(mode==1){while(k<num){fun(a,b,x,y,&max,N);k++;if(k<10)printf("%d ",k);elseprintf("%d ",k);for(i=1;i<=N;i++)printf("x[%d]=%-10.7lf",i,y[i]);printf("\n");}}else{while(max>e){fun(a,b,x,y,&max,N);k++;if(k<10)printf("%d ",k);elseprintf("%d ",k);for(i=1;i<=N;i++)printf("x[%d]=%-10.7lf",i,y[i]);printf("\n");}}getchar();getchar();}共轭梯度法c语言源程序:#include<stdio.h>#include<math.h>#define N 100void fun1(double a[][N],double *b,double *c){int i,j;double sum=0;for(i=1;i<N;i++){sum=0;for(j=1;j<N;j++)sum+=a[i][j]*b[j];c[i]=sum;}}double fun2(double *c,double *d){int i;double sum=0;for(i=1;i<N;i++)sum+=c[i]*d[i];return sum;}void main(){int i,j,mode,n,num,k=0;double e,a[N][N]={0};double b[N]={0},x0[N]={0},x[N]={0},r[N]={0},d[N]={0};double Beta,Lambda,max=0,c1[N]={0},c2[N]={0};printf("欢迎来到共轭梯度法解线性方程组，请先确保该法收敛，按提示操作*^_^*\n");printf("请输入结束条件，迭代次数限制请按1，误差限制请按2:");scanf("%d",&mode);if(mode==2){printf("请输入误差限e:");scanf("%lf",&e);}else{printf("请输入迭代次数：");scanf("%d",&num);}printf("请输入方程阶数：");scanf("%d",&n);printf("请输入系数矩阵：\n");for(i=1;i<=n;i++)for(j=1;j<=n;j++)scanf("%lf",&a[i][j]);printf("请输入向量b:\n");for(i=1;i<=n;i++)scanf("%lf",&b[i]);printf("请输入初始解x:");for(j=1;j<=n;j++){scanf("%lf",&x[j]);x0[j]=x[j];}fun1(a,x,c1);for(i=1;i<=n;i++){d[i]=b[i]-c1[i];r[i]=d[i];max=r[1];if(fabs(r[i])>max)max=fabs(r[i]);}fun1(a,d,c1);Lambda=fun2(r,d)/fun2(d,c1);printf("%d ",k);for(i=1;i<=n;i++)printf("x[%d]=%-10.7lf",i,x0[i]); printf("\n");k++;printf("%d ",k);for(i=1;i<=n;i++){x[i]=x[i]+Lambda*d[i];printf("x[%d]=%-10.7lf",i,x[i]);}printf("\n");if(mode==1){while(k<num){k++;fun1(a,x,c1);for(i=1;i<=n;i++){r[i]=b[i]-c1[i];max=r[1];if(fabs(r[i])>max)max=fabs(r[i]);}fun1(a,d,c2);Beta=-fun2(r,c2)/fun2(d,c2);for(i=1;i<=n;i++)d[i]=r[i]+Beta*d[i];fun1(a,d,c2);Lambda=fun2(r,d)/fun2(d,c2);if(k<10)printf("%d ",k);elseprintf("%d ",k);for(i=1;i<=n;i++){x[i]+=Lambda*d[i];printf("x[%d]=%-10.7lf",i,x[i]);}printf("\n");}}else{while(max>e){k++;fun1(a,x,c1);for(i=1;i<=n;i++){r[i]=b[i]-c1[i];max=r[1];if(fabs(r[i])>max)max=fabs(r[i]);}fun1(a,d,c2);Beta=-fun2(r,c2)/fun2(d,c2);for(i=1;i<=n;i++)d[i]=r[i]+Beta*d[i];fun1(a,d,c2);Lambda=fun2(r,d)/fun2(d,c2);if(k<10)printf("%d ",k);elseprintf("%d ",k);for(i=1;i<=n;i++){x[i]+=Lambda*d[i];printf("x[%d]=%-10.7lf",i,x[i]);}printf("\n");}}getchar();getchar();}程序代码（matlab）（1）Jacobi迭代法MatLab源程序。