计算方法上机作业
（完整版）数值计算方法上机实习题答案
（完整版）数值计算⽅法上机实习题答案1．设?+=105dx xx I nn ，（1）由递推公式nI I n n 151+-=-，从0I 的⼏个近似值出发，计算20I ；解：易得：0I =ln6-ln5=0.1823, 程序为：I=0.182; for n=1:20I=(-5)*I+1/n; end I输出结果为：20I = -3.0666e+010 （2）粗糙估计20I ，⽤nI I n n 515111+-=--，计算0I ；因为 0095.056 0079.01020201020≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2120=+=I 程序为：I=0.0087; for n=1:20I=(-1/5)*I+1/(5*n); end I0I = 0.0083（3）分析结果的可靠性及产⽣此现象的原因(重点分析原因)。
⾸先分析两种递推式的误差；设第⼀递推式中开始时的误差为000I I E '-=，递推过程的舍⼊误差不计。
并记nn n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。
因为=20E 20020)5(I E >>-，所此递推式不可靠。
⽽在第⼆种递推式中n n E E E )51(5110-==-=Λ，误差在缩⼩，所以此递推式是可靠的。
出现以上运⾏结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。
2．求⽅程0210=-+x e x的近似根，要求41105-+?<-k k x x ，并⽐较计算量。
（1）在[0，1]上⽤⼆分法；程序：a=0;b=1.0;while abs(b-a)>5*1e-4 c=(b+a)/2;if exp(c)+10*c-2>0 b=c; else a=c; end end c结果：c =0.0903（2）取初值00=x ，并⽤迭代1021x k e x -=+；程序：x=0; a=1;while abs(x-a)>5*1e-4 a=x;x=(2-exp(x))/10; end x结果：x =0.0905（3）加速迭代的结果；程序：x=0; a=0;b=1;while abs(b-a)>5*1e-4 a=x;y=exp(x)+10*x-2; z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x); b=x; end x结果：x =0.0995（4）取初值00=x ，并⽤⽜顿迭代法；程序：x=0; a=0;b=1;while abs(b-a)>5*1e-4 a=x;x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x; end x结果： x =0.0905（5）分析绝对误差。
计算方法上机题答案
2.用下列方法求方程e^x+10x-2=0的近似根，要求误差不超过5*10的负4次方，并比较计算量（1）二分法（局部，大图不太看得清，故后面两小题都用局部截图）（2）迭代法（3）牛顿法顺序消元法#include<stdio.h>#include<stdlib.h>#include<math.h>int main(){ int N=4,i,j,p,q,k; double m;double a[4][5]; double x1,x2,x3,x4; for (i=0;i<N ;i++ )for (j=0;j<N+1; j++ )scanf("%lf",&a[i][j]);for(p=0;p<N-1;p++) {for(k=p+1;k<N;k++){m=a[k][p]/a[p][p];for(q=p;q<N+1;q++)a[k][q]=a[k][q]-m*a[p][q];}}x4=a[3][4]/a[3][3];x3=(a[2][4]-x4*a[2][3])/a[2][2];x2=(a[1][4]-x4*a[1][3]-x3*a[1][2])/a[1][1];x1=(a[0][4]-x4*a[0][3]-x3*a[0][2]-x2*a[0][1])/a[0][0];printf("%f,%f,%f,%f",x1,x2,x3,x4);scanf("%lf",&a[i][j]); （这一步只是为了看到运行的结果）}运行结果列主元消元法function[x,det,flag]=Gauss(A,b)[n,m]=size(A);nb=length(b);flag='OK';det=1;x=zeros(n,1);for k=1:n-1 max1=0;for i=k:nif abs(A(i,k))>max1max1=abs(A(i,k));r=i;endendif max1<1e-10flag='failure';return;endif r>kfor j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det; endfor i=k+1:nm=A(i,k)/A(k,k);for j=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A(n,n)if abs(A(n,n))<1e-10flag='failure';return;endx(n)=b(n)/A(n,n);for k=n-1:-1：1for j=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);end运行结果：雅可比迭代法function y=jacobi(a,b,x0) D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;while norm(y-x0)>1e-4 x0=y;y=B*x0+f;n=n+1;endyn高斯赛德尔迭代法function y=seidel(a,b,x0) D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1;while norm(y-x0)>10^(-4) x0=y;y=G*x0+f;n=n+1;endynSOR迭代法function y=sor(a,b,w,x0)D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);lw=(D-w*L)\((1-w)*D+w*U); f=(D-w*L)\b*w;y=lw*x0+f;n=1;while norm(y-x0)>10^(-4) x0=y;y=lw*x0+f;n=n+1;endyn1.分段线性插值：function y=fdxx(x0,y0,x)p=length(y0);n=length(x0);m=length(x);for i=1:mz=x(i);for j=1:n-1if z<x0(j+1)break;endendy(i)= y0(j)*(z-x0(j+1))/(x0(j)-x0(j+1))+y0(j+1)*(z-x0(j))/(x0(j+1)-x0(j));fprintf('y(%d)=%f\nx1=%.3fy1=%.3f\nx2=%.3fy2=%.3f\n\n',i,y(i),x0(j),y0(j),x0(j+1),y0(j+1));endend结果0.39404 0.38007 0.356932.分段二次插值：function y=fdec(x0,y0,x)p=length(y0);n=length(x0);m=length(x);for i=1:mz=x(i);for j=1:n-1if z<x0(j+1)break;endendif j<2j=j+1;elseif (j<n-1)if (abs(x0(j)-z)>abs(x0(j+1)-z))j=j+1;elseif ((abs(x0(j)-z)==abs(x0(j+1)-z))&&(abs(x0(j-1)-z)>abs(x0(j+2)-z)))j=j+1;endendans=0.0;for t=j-1:j+1a=1.0;for k=j-1:j+1if t~=ka=a*(z-x0(k))/(x0(t)-x0(k));endendans=ans+a*y0(t);endy(i)=ans;fprintf('y(%d)=%f\n x1=%.3f y1=%.3f\n x2=%.3f y2=%.3f\n x3=%.3f y3=%.3f\n\n',i,y(i),x0(j-1),y0(j-1),x0(j),y0(j),x0(j+1),y0(j+1));endend结果为0.39447 0.38022 0.357253.拉格朗日全区间插值function y=lglr(x0,y0,x)p=length(y0);n=length(x0);m=length(x);for t=1:mans=0.0;z=x(t);for k=1:np=1.0;for q=1:nif q~=kp=p*(z-x0(q))/(x0(k)-x0(q));endendans=ans+p*y0(k);endy(t)=ans;fprintf('y(%d)=%f\n',t,y(t));endend结果为0.39447 0.38022 0.35722function [p,S,mu] = polyfit(x,y,n)if ~isequal(size(x),size(y))errorendx = x(:);y = y(:);if nargout > 2mu = [mean(x); std(x)];x = (x - mu(1))/mu(2);endV(:,n+1) = ones(length(x),1,class(x));for j = n:-1:1V(:,j) = x.*V(:,j+1);end[Q,R] = qr(V,0);ws = warning;p = R\(Q'*y); warning(ws);if size(R,2) > size(R,1)warning;elseif warnIfLargeConditionNumber(R)if nargout > 2warning;elsewarning;endendif nargout > 1r = y - V*p;S.R = R;S.df = max(0,length(y) - (n+1));S.normr = norm(r);endp = p.';function flag = warnIfLargeConditionNumber(R) if isa(R, 'double')flag = (condest(R) > 1e+10);elseflag = (condest(R) > 1e+05);endx=[1,3,4,5,6,7,8,9,10];y=[10,5,4,2,1,1,2,3,4];p=polyfit(x,y,2);y=poly2sym(p);a=-p(1)/p(2)*0.5;xa=-p(2)/(2*p(1));min=(4*p(1)*p(3)-p(2)^2)/(4*p(1)); yxamin运行截图运行结果function [I] = CombineTraprl(f,a,b,eps)if(nargin ==3)eps=1.0e-4;endn=1;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>epsn=n+1;h=(b-a)/n;I1=I2;I2=0;for i=0:n-1x=a+h*i;x1=x+h;I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym( f),findsym(sym(f)),x1));endendI=I2;程序：>> [q]=CombineTraprl('(1-exp(-x))^0.5/x'10^-12,1)结果：q =1.8521欧拉方法function [x,y]=euler(fun,x0,xfinal,y0,n);if nargin<5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;for i=1:nx(i+1)=x(i)+h;y(i+1)=y(i)+h*feval(fun,x(i),y(i));end程序：[x,y]=euler('doty',0,1,1,10)结果：改进欧拉方法function [x,y]=eulerpro(fun,x0,xfinal,y0,n); if nargin<5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;x(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i),y(i));y2=y(i)+h*feval(fun,x(i+1),y1);y(i+1)=(y1+y2)/2;end程序：>> [x,y]=eulerpro('doty',0,1,1,10)结果：经典RK法function [x,y]=RungKutta4(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1k1=feval(dyfun,x(n),y(n));k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1);k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2);k4=feval(dyfun,x(n+1),y(n)+h*k3);y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6; end程序：dyfun=inline('(2*x*y^-2)/3');[x,y]=RungKutta4(dyfun,[0,1],1,0.1)标准函数x(1)=0;h=0.1;for i=1:10y(i)=(1+x(i)^2)^(1/3); endxY结果：。
计算方法上机作业——龙格库塔法
计算方法上机作业——龙格库塔法龙格库塔法（Runge-Kutta method）是一种常用于求解常微分方程（Ordinary Differential Equation，ODE）的数值解法。
它是由德国数学家卡尔·龙格（Carl Runge）和马丁·威尔海姆·库塔（Martin Wilhelm Kutta）在20世纪初提出的。
龙格库塔法的基本思想是通过数值逼近来计算微分方程的近似解。
在讲解龙格库塔法之前，我们先来简单回顾一下ODE的一阶常微分方程的基本形式：y′(y)=y(y,y)其中，y(y,y)是已知函数。
龙格库塔法的核心是使用差分逼近计算函数的斜率。
假设我们要求解的方程为：y′(y)=y(y,y),y(y)=y₀所需计算的点为y₀,y₁,...,yy，对应的函数值为y₀,y₁,...,yy，其中y是步长的个数。
龙格库塔法通过递推关系式来计算估计值，并不断更新当前点的函数值。
接下来以龙格库塔法的经典四阶形式为例进行说明。
该方法的基本方程如下：yy+1=yy+(y₁+2y₂+2y₃+y₄)/6y₁=ℎy(yy,yy)y₂=ℎy(yy+ℎ/2,yy+y₁/2)y₃=ℎy(yy+ℎ/2,yy+y₂/2)y₄=ℎy(yy+ℎ,yy+y₃)其中y表示当前步骤，ℎ表示步长，yy表示当前点的函数值，y₁,y₂,y₃和y₄则表示对应的斜率。
使用龙格库塔法，我们可以通过不断递归计算来求得指定区间（例如[y,y]）上的函数值。
具体步骤如下：1.确定求解区间[y,y]和初始点(y₀,y₀)以及步长ℎ。
2.初始化：设置yy=y₀，yy=y₀。
3.对所有y=0,...,y−1：计算y₁,y₂,y₃和y₄，根据上述递推关系式。
根据递推关系式计算yy+1更新当前点的函数值，即yy+1=y(yy+1)。
更新当前点的y值，即yy+1=yy+ℎ。
4.返回结果：最终求得的函数值。
需要注意的是，选择适当的步长对最终结果的精度和计算效率都有重要影响。
数值计算方法上机实验报告
数值计算方法上机实验报告
一、实验目的
本次实验的主要目的是熟悉和掌握数值计算方法，学习梯度下降法的
原理和实际应用，熟悉Python语言的编程基础知识，掌握Python语言的
基本语法。
二、设计思路
本次实验主要使用的python语言，利用python下的numpy，matplotlib这两个工具，来实现数值计算和可视化的任务。
1. 首先了解numpy的基本使用方法，学习numpy的矩阵操作，以及numpy提供的常见算法，如矩阵分解、特征值分解等。
2. 在了解numpy的基本操作后，可以学习matplotlib库中的可视化
技术，掌握如何将生成的数据以图表的形式展示出来。
3. 接下来就是要学习梯度下降法，首先了解梯度下降法的主要原理，以及具体的实际应用，用python实现梯度下降法给出的算法框架，最终
可以达到所期望的优化结果。
三、实验步骤
1. 熟悉Python语言的基本语法。
首先是熟悉Python语言的基本语法，学习如何使用Python实现变量
定义，控制语句，函数定义，类使用，以及面向对象编程的基本概念。
2. 学习numpy库的使用方法。
其次是学习numpy库的使用方法，学习如何使用numpy库构建矩阵，学习numpy库的向量，矩阵操作，以及numpy库提供的常见算法，如矩阵分解，特征值分解等。
3. 学习matplotlib库的使用方法。
东南大学计算方法与实习上机实验一
东南大学计算方法与实习实验报告学院：电子科学与工程学院学号：06A*****姓名：***指导老师：***实习题14、设S N=Σ(1)编制按从大到小的顺序计算S N的程序；(2)编制按从小到大的顺序计算S N的程序；(3)按两种顺序分别计算S1000，S10000，S30000，并指出有效位数。
解析：从大到小时，将S N分解成S N-1=S N-,在计算时根据想要得到的值取合适的最大的值作为首项；同理从小到大时，将S N=S N-1+ ，则取S2=1/3。
则所得式子即为该算法的通项公式。
(1)从大到小算法的C++程序如下：/*从大到小的算法*/#include<iostream>#include<iomanip>#include<cmath>using namespace std;const int max=34000; //根据第(3)问的问题，我选择了最大数为34000作为初值void main(){int num;char jus;double cor,sub;A: cout<<"请输入你想计算的值"<<'\t';cin>>num;double smax=1.0/2.0*(3.0/2.0-1.0/max-1.0/(max+1)),temps;double S[max];// cout<<"s["<<max<<"]="<<setprecision(20)<<smax<<'\n';for(int n=max;n>num;){temps=smax;S[n]=temps;n--;smax=smax-1.0/((n+1)*(n+1)-1.0);}cor=1.0/2.0*(3.0/2.0-1.0/num-1.0/(num+1.0)); //利用已知精确值公式计算精确值sub=fabs(cor-smax); //double型取误差的绝对值cout<<"用递推公式算出来的s["<<n<<"]="<<setprecision(20)<<smax<<'\n';cout<<"实际精确值为"<<setprecision(20)<<cor<<'\n';cout<<"则误差为"<<setprecision(20)<<sub<<'\n';cout<<"是否继续计算S[N],是请输入Y，否则输入N！"<<endl;cin>>jus;if ((int)jus==89||(int)jus==121) goto A;}(2)从小到大算法的C++程序如下：/*从小到大的算法*/#include<iostream>#include<iomanip>#include<cmath>using namespace std;void main(){int max;A: cout<<"请输入你想计算的数，注意不要小于2"<<'\t';cin>>max;double s2=1.0/3.0,temps,cor,sub;char jus;double S[100000];for(int j=2;j<max;){temps=s2;S[j]=temps;j++;s2+=1.0/(j*j-1.0);}cor=1.0/2.0*(3.0/2.0-1.0/j-1.0/(j+1.0)); //利用已知精确值公式计算精确值sub=fabs(cor-s2); //double型取误差的绝对值cout<<"用递推公式算出来的s["<<j<<"]="<<setprecision(20)<<s2<<'\n';cout<<"实际精确值为"<<setprecision(20)<<cor<<'\n';cout<<"则误差为"<<setprecision(20)<<sub<<'\n';cout<<"是否继续计算S[N],是请输入Y，否则输入N！"<<endl;cin>>jus;if ((int)jus==89||(int)jus==121) goto A;}(3)(注：因为程序中setprecision(20)表示输出数值小数位数20，则程序运行时所得到的有效数字在17位左右)ii.选择从小到大的顺序计算S1000、S10000、S30000的值需要计算的项S1000S10000S30000计算值0.74900049950049996 0.74966672220370571 0.74996666722220728实际精确值0.74900049950049952 0.74990000499950005 0.74996666722220373误差 4.4408920985006262*10-16 5.6621374255882984*10-15 3.5527136788005009*10-15有效数字17 17 17附上部分程序运行图：iii.实验分析通过C++程序进行计算验证采用从大到小或者从小到大的递推公式算法得到的数值基本稳定且误差不大。
东南大学计算方法上机报告实验报告完整版
实习题11. 用两种不同的顺序计算644834.11000012≈∑=-n n，试分析其误差的变化解：从n=1开始累加，n 逐步增大，直到n=10000；从n=10000开始累加，n 逐步减小，直至1。
算法1的C 语言程序如下： #include<stdio.h> #include<math.h> void main() { float n=0.0; int i; for(i=1;i<=10000;i++) { n=n+1.0/(i*i); } printf("%-100f",n); printf("\n"); float m=0.0; int j; for(j=10000;j>=1;j--) { m=m+1.0/(j*j); } printf("%-7f",m); printf("\n"); }运行后结果如下：结论： 4.设∑=-=Nj N j S 2211，已知其精确值为)11123(21+--N N 。
1）编制按从大到小的顺序计算N S 的程序； 2）编制按从小到大的顺序计算N S 的程序；3）按2种顺序分别计算30000100001000,,S S S ，并指出有效位数。
解：1）从大到小的C语言算法如下：#include<stdio.h>#include<math.h>#include<iostream>using namespace std;void main(){float n=0.0;int i;int N;cout<<"Please input N"<<endl;cin>>N;for(i=N;i>1;i--){n=n+1.0/(i*i-1);N=N-1;}printf("%-100f",n);printf("\n");}执行后结果为：N=2时，运行结果为：N=3时，运行结果为：N=100时，运行结果为：N=4000时，运行结果为：2）从小到大的C语言算法如下：#include<stdio.h>#include<math.h>#include<iostream>using namespace std;void main(){float n=0.0;int i;int N;cout<<"Please input N"<<endl;cin>>N;for(i=2;i<=N;i++){n=n+1.0/(i*i-1);}printf("%-100f",n);printf("\n");}执行后结果为：N=2时，运行结果为：N=3时，运行结果为：N=100时，运行结果为：N=4000时，运行结果为：结论：通过比较可知：N 的值较小时两种算法的运算结果相差不大，但随着N 的逐渐增大，两种算法的运行结果相差越来越大。
应用计算方法
结果分析：比较发现，经过两种改进迭代法，求重根时迭代速度明显加快。
3-4 试验目的体验 Steffensen’s method 加速技巧 试验内容：先用 Newton 法求解方程 x-tanx=0 再用 Steffensen 法求解，比较迭代步数。精确到 10-4。 迭代公式： P(k+1)=P(k)-(P(k)-tan(P(k)))/(1-(sec(P（k）))^2) 初值 P0=1 Newton 法： function [ p,k ]=fnewton( p0,max,tol ) for k=1:max p=p0-(p0-tan(p0))/(1-(sec(p0))^2); if abs(p-p0)<tol break; end p0=p; end disp(p); disp(k) % fnewton( 1,100,10^(-4) ) Steffensen 法： function [ p1,k ]=steffensen( p0,max,tol ) n=1; p(1)=p0; while n<=max for k=1:2 p(k+1)=p(k)-(p(k)-tan(p(k)))/(1-(sec(p(k)))^2); end p1=p(1)-(p(2)-p(1))^2/(p(3)-2*p(2)+p(1)); f0=p1-(p1-tan(p1))/(1-(sec(p1))^2); if abs(f0)<tol break; end n=n+1; p(1)=p1; end disp(p1); disp(n) % steffensen( 1,100,10^(-4) ) >> steffensen( 1,100,10^(-4) ) -1.3387e-07 3
3-2 试验目的:考察 Newton 法求单根的收敛速度 应用 Newton 迭代法求 3-1 中的方程，并与 3-1 中的迭代法相比较，考察收敛速度。精确到时 10-4 迭代公式： P(k+1)= P(k)-(2P(k)-eP(k)+3)/(2- eP(k)) 初值 P0=1 和-1。 Newton 法: function [ p,k ] = fnewton( p0,max,tol ) for k=1:max p=p0-(2*p0-exp(p0)+3)/(2-exp(p0)); if abs(p-p0)<tol break; end
计算方法上机作业
计算方法上机报告姓名：学号：班级：上课班级：说明：本次上机实验使用的编程语言是Matlab 语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。
1. 对以下和式计算：∑∞⎪⎭⎫ ⎝⎛+-+-+-+=0681581482184161n n n n S n，要求：①若只需保留11个有效数字，该如何进行计算； ②若要保留30个有效数字，则又将如何进行计算； （1） 算法思想1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为：142111416818485861681n n na n n n n n ε⎛⎫=---<< ⎪+++++⎝⎭；2、为了保证计算结果的准确性，写程序时，从后向前计算；3、使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数)（2）算法结构1.;0=s⎪⎭⎫⎝⎛+-+-+-+=681581482184161n n n n t n ；2.for 0,1,2,,n i =⋅⋅⋅if 10mt -≤end;3.for ,1,2,,0n i i i =--⋅⋅⋅;s s t =+（3）Matlab 源程序clear; %清除工作空间变量 clc; %清除命令窗口命令m=input('请输入有效数字的位数m='); %输入有效数字的位数s=0;for n=0:50t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));if t<=10^(-m) %判断通项与精度的关系break;endend;fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值for i=n-1:-1:0t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));s=s+t; %求和运算ends=vpa(s,m) %控制s的精度（4）结果与分析当保留11位有效数字时，需要将n值加到n=7，s =3.1415926536；当保留30位有效数字时，需要将n值加到n=22,s =3.14159265358979323846264338328。
计算方法上机实习作业
如图所示，在两端出现龙格现象，并且阶数越高龙格现象越明显。 若在区间［0，5］进行插值，是否还存在龙格现象？
Matlab 程序如下： x0=0:0.1:5; y0=1./(1+x0.^2); plot(x0,y0,'r'); hold on; for i=3:2:11 x=linspace(-0,5,i); y=1./(1+x.^2); p=polyfit(x,y,i-1); xx=0:0.1:5; yy=polyval(p,xx); plot(xx,yy,'b'); hold on; grid on; end;
0
1
2
3
4
Matlab 程序如下： clear;clc; x=0:0.1:2.5; y=sin(x); y1=x; y2=x-x.^3/6; y3=x-x.^3/6+x.^5/120; plot(x,y,'b',x,y1,'g',x,y2,'r',x,y3,'k')
输出结果如图：
1.500000000000000 1.414213562374690 1.414213562373095
1.416666666666667 1.414213562373095
2.000000000000000 1.732050810014727 1.732050807568877
1.750000000000000 1.732050807568877
x 2 n 1 ，编写程序实 ( 2n 1)! k 1 现如下功能：对 n=1，2，3，绘制出区间[0，2.5 ]上的近似函数的图形以及 sin x 本身的图
4.由正弦函数的 Tylor 展开式取前 n 项，得 sin x ( 1) n1
计算方法 上机操作(C语言)
#include "stdlib.h"#include "math.h"#include "stdio.h"#define N 4intagaus(double a[N][N],double b[N],int n);//高斯全主元法求deltaX，具体见程序中注释void GetFun(double x[N],double b[N]);//F(X),结果放在b里void GetJacobi(double x[N],double J[N][N]);//计算F'(X)，结果放在J里double GetDeltaX(double J[N][N],double b[N],double deltaX[N]);//求deltaX主程序，调用agaus voidUpdateX(double x[N],double deltaX[N]);//xk+1=xk+deltaXvoid SolveNonlinear(double x[N],double eps,int k);//非线性方程组主程序,迭代终止条件是最大迭代次数或者绝对值范数（街区距离）main(){inti,k;doubleeps;double x[4]={1,1.0,1.0,1.0};double b[N];eps=0.0000001; k=100;SolveNonlinear(x,eps,k);GetFun(x,b);printf("\n");printf("\n");for (i=0; i<N; i++)printf("x(%d)=%13.7e,F(x)=%13.7e\n",i,x[i],b[i]);printf("\n");getchar();return 0;}voidSolveNonlinear(double x0[N],double eps,int k){intnIter=0;double b[N],J[N][N],deltaX[N],NormDeltaX;while(1){nIter++;GetFun(x0,b);GetJacobi(x0,J);NormDeltaX=GetDeltaX(J,b,deltaX);UpdateX(x0,deltaX);if ((nIter>k)||(NormDeltaX<eps))break;}}void GetFun(x,y)//F(X),结果放在b里double x[],y[];{y[0]=x[0]*x[0]+x[1]*x[1]+x[2]*x[2]+x[3]*x[3]-1.0;y[1]=2.0*x[0]*x[0]+x[1]*x[1]-4.0*x[2]+2*x[3]*x[3];y[2]=3.0*x[0]*x[0]-4.0*x[1]+x[2]*x[2]+3*x[3]*x[3];y[3]=4.0*x[0]+x[1]*x[1]+2*x[2]*x[2]+5*x[3]*x[3]-3.0;return;}void GetJacobi(x,J)//计算F'(X)，结果放在J里double x[N],J[N][N];{J[0][0]=2*x[0];J[0][1]=2*x[1];J[0][2]=2*x[2];J[0][3]=2*x[3];J[1][0]=4*x[0];J[1][1]=2*x[1];J[1][2]=-4;J[1][3]=4*x[3];J[2][0]=6*x[0];J[2][1]=-4;J[2][2]=2*x[2];J[2][3]=6*x[3];J[3][0]=4;J[3][1]=2*x[1];J[3][2]=4*x[2];J[3][3]=10*x[3];return;}intagaus(a,b,n)//高斯全主元法解方程，结果放在b里int n;double a[N][N],b[N];{int *js,l,k,j,i,is;double d,t;js=malloc(n*sizeof(int));//记载列交换标号l=1;for (k=0;k<=n-2;k++){d=0.0;for (i=k;i<=n-1;i++)//求全主元{for (j=k;j<=n-1;j++){t=fabs(a[i][j]);if (t>d){d=t;js[k]=j;//记住列is=i;//记住行}}}if (d+1.0==1.0)l=0;else{/* if (js[k]!=k)//列交换{for (i=0;i<=n-1;i++){t=a[i][k]; a[i][k]=a[i][js[k]]; a[i][js[k]]=t;}}*/if (is!=k){for (j=k;j<=n-1;j++)//行交换{t=a[k][j]; a[k][j]=a[is][j]; a[is][j]=t;}t=b[k]; b[k]=b[is]; b[is]=t;//b交换}}if (l==0)//最大主元为零不能计算出正确结果{free(js);printf("fail\n");return(0);}//以下为消元过程d=a[k][k];a[k][k]=1;//这一步是否要无关紧要，只要在回代时不除以该项即可。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
方程求根一、目的和意义非线性方程在科学研究与工程实践中广泛出现，例如，优化问题、特征值问题、微分方程问题等。
但是，除少量方程外，大多数非线性方程求根相当困难，常见的几个简单、有效的数值求根方法，包括二分法、迭代法、牛顿法、割线法等，本实验旨在比较各种算法的计算性能和使用范围。
二、计算公式1.二分法2.不动点迭代法三、结构程序设计代码1.二分法1).定义所求解函数function [ y ] = f( x )y = x^3 + 4*x^2 - 10;end2).执行算法%初始化，设置区间端点a、b，误差限tola = 1;b = 2; tol = 0.5*10^(-6);k = 0; fa = f(a);%设置最大二分次数为30for k = 1:50p = (a + b)/2; fp = f(p);if(fp == 0 || (b - a)/2 < tol)breakendif(fa * fp < 0)b = p;elsea = p;enddisp('近似解p = ');disp(vpa(p,10)); disp('迭代次数k = ');disp(k);end2.不动点迭代法1).定义不动点方程g(x)function [ y ] = g( x )y = x^3 + 4*x^2 + x - 10;end2).执行算法%初始化，设置误差限，设置初值p0tol = 0.5*10^(-6);k = 1; p0 = 1.5;%迭代次数为10次while k <= 10p = g(p0);if abs(p - p0) < tolbreakenddisp('近似解p = ');disp(vpa(p,10)); disp('迭代次数k = ');disp(k);k = k + 1; p0 = p;end四、结果及其讨论1.二分法结果由于结果较长，只取了一部分，从图中可以看出，迭代20次可得到误差限范围内的近似解p=1.36522。
2.不动点迭代法结果由结果可以看出，函数不收敛，当迭代次数为6时，已经近似趋向于无穷大。
3.讨论上面实验通过两种方法近似计算了函数的解，二分法可在迭代20次的情况下求出近似解，而由于所设的不动点方程发散，所以使用不动点迭代法无法求出近似解。
因此，在实际计算过程中，使用不动点迭代法时，在建立不动点方程后，要分析其收敛性。
线性方程组的直接解法一、目的和意义在科学和工程技术中，很多问题可归结为求解大型线性方程组。
很多非线性问题也可通过局部线性化而转化为线性问题。
因此，线性问题的求解构成了解决各类计算问题的基础。
线性方程组的数值解法分为直接解法和间接解法两大类。
直接解法是指通过有限步运算后，将方程组加工为某个三角方程组或者对角方程组。
直接法常用于求解中小型的线性方程组及某些特殊的大型方程组。
二、计算公式1.LU分解法(doolittle分解)三、结构程序设计代码1.LU分解法题1%初始化，矩阵A，4阶，列向量b，L为单位对角矩阵A = [4 2 1 5;8 7 2 10;4 8 3 6;12 6 11 20]; n=4; b = [-2 -7 -7 -3]';L = eye(n); U(1,:) = A(1,:);for k = 2:n%判断U矩阵的对角线上是否有0，如果有0，则不能分解if U(k-1,k-1) == 0disp('分解失败');returnend%具体分解过程L(k:n, k-1) = A(k:n, k-1)/U(k-1, k-1);U(k, k:n) = A(k, k:n) - L(k, 1:k-1) * U(1:k-1, k:n);if k < nA(k+1:n, k) = A(k+1:n, k) - L(k+1:n, 1:k-1) * U(1:k-1,k);endend%利用Ly=b，算出y向量y = (L^-1) * b;%利用Ux=y，算出解向量xx = (U^-1) * y题2-1%初始化，矩阵A，n阶，列向量b，L为单位对角矩阵n=5; A = zeros(n,n); b = zeros(n,1);%循环构建矩阵A和向量bfor i = 1:nA(i,i) = 2;if i < nA(i,i+1) = 1; A(i+1,i) = 1;endif i == 1b(i,1) = -7;elseb(i,1) = -5;endendL = eye(n); U(1,:) = A(1,:);for k = 2:n%判断U矩阵的对角线上是否有0，如果有0，则不能分解if U(k-1,k-1) == 0disp('分解失败');returnend%具体分解过程L(k:n, k-1) = A(k:n, k-1)/U(k-1, k-1);U(k, k:n) = A(k, k:n) - L(k, 1:k-1) * U(1:k-1, k:n);if k < nA(k+1:n, k) = A(k+1:n, k) - L(k+1:n, 1:k-1) * U(1:k-1,k);endendLU%利用Ly=b，算出y向量y = (L^-1) * b;%利用Ux=y，算出解向量xx = (U^-1) * y题2-2、2-3这两个程序与2-1的程序几乎一模一样，只是将n的值变为10和100，再进行计算，为了节省篇幅，这里就不再列出程序了。
四、结果及其讨论1.LU分解法结果题1题2-1题2-2题2-3为了节省篇幅，这里不再列出L和U矩阵，只将解列出。
x =-4.45541.9109-4.36631.8218-4.27721.7327-4.18811.6436-4.09901.5545 -4.0099 1.4653 -3.9208 1.3762 -3.8317 1.2871 -3.7426 1.1980 -3.6535 1.1089 -3.5644 1.0198 -3.4752 0.9307 -3.3861 0.8416 -3.2970 0.7525 -3.2079 0.6634 -3.1188 0.5743 -3.0297 0.4851 -2.9406 0.3960 -2.8515 0.3069 -2.7624 0.2178 -2.6733 0.1287 -2.5842 0.0396 -2.4950 -0.0495 -2.4059 -0.1386 -2.3168 -0.2277 -2.2277 -0.3168 -2.1386-2.0495 -0.4950 -1.9604 -0.5842 -1.8713 -0.6733 -1.7822 -0.7624 -1.6931 -0.8515 -1.6040 -0.9406 -1.5149 -1.0297 -1.4257 -1.1188 -1.3366 -1.2079 -1.2475 -1.2970 -1.1584 -1.3861 -1.0693 -1.4752 -0.9802 -1.5644 -0.8911 -1.6535 -0.8020 -1.7426 -0.7129 -1.8317 -0.6238 -1.9208 -0.5347 -2.0099 -0.4455 -2.0990 -0.3564 -2.1881 -0.2673 -2.2772 -0.1782-0.0891-2.45543.讨论LU分解在本质上是高斯消元法的一种表达形式。
实质上是将A通过初等行变换变成一个上三角矩阵，其变换矩阵就是一个单位下三角矩阵。
这正是所谓的杜尔里特算法（Doolittle algorithm）：从下至上地对矩阵A做初等行变换，将对角线左下方的元素变成零，然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积的逆就是L矩阵，它也是一个单位下三角矩阵。
这类算法的复杂度一般在(三分之二的n三次方) 左右。
在计算第2题的过程中，当把n=5的情况算完之后，将n改为10继续进行计算，发现matlab会报错。
经过反复思考发现，因为程序只修改了n的值，其他参数（如矩阵）的名称并未进行修改，matlab也会将这些矩阵的值存起来，所以在计算时候会报错。
因此，当运行完一次计算会，应该在matlab主界面的命令行输入clear指令，把前面计算的值都清楚掉，这样就不会影响到之后的计算了。
线性方程组的迭代解法一、目的和意义直接法对求解中小型方程组很有效，但当线性方程组的变量个数较多时，直接解法由于存储量和计算量太大就不是用来。
此外，直接法也不适用于求解非线性方程组。
求解线性方程组的迭代法的基本思想是将“复杂”化归“简单”的重复。
具体来讲，迭代法是一种逐次逼近的方法，从线性方程组的一个假设解向量开始，通过来自线性方程组的迭代格式，产生一系列近似解向量，最后求出满足精度要求的近似解。
二、计算公式SOR方法三、结构程序设计代码SOR方法%初始化，矩阵A，列向量b，初始向量X0，松弛因子w，最大迭代次数NA = [10 -1 -2;-1 10 -2;-1 -1 5]; b = [7.2 8.3 4.2]';X0 = [1 1 1]'; w = 0.5; N = 30;X = X0;for K = 1:N%使i=1:3，循环得出各次迭代的X(1)，X(2)，X(3)的值，即得出各次迭代后的X向量for i = 1:3X(i) = w * (b(i) - A(i,:)*X) / A(i,i) + X(i);endX0 = X;%每次循环过后，打印迭代的次数和解向量X的值disp(K);disp(X);end另外，为了验证不同松弛因此对SOR迭代法的影响，此后我们将松弛因子w设为1，即Gauss-Seidel迭代法，再进行计算，比较两次的结果。
四、结果及其讨论1.SOR迭代法，松弛因子w=0.5由结果可以看出，计算迭代到第20次的时候，就可以得到较为精确的结果，解向量为X=[1.1 1.2 1.3]T。
2.SOR迭代法，松弛因子w=1，即Gauss-Seidel迭代法由结果可以看出，计算迭代到第6次的时候，就可以得到较为精确的结果。
3.讨论通过实验可知，对于SOR迭代法，选取不同的松弛因子，在计算的过程中迭代次数会有不同。
因此，应用SOR法的关键是选择合适的松弛因子。
但在一般情况下，如何选取最佳的松弛因子，以提高迭代向量序列的收敛速度是比较复杂的。
通常是在(0,2)中选取几个数据，比较迭代向量序列的收敛速度，选择一个最快的数据就可以了，当然这种方法只是一种试探性的方法。