数值计算方法上机实验报告
西南交通大学数值分析上机实验报告
数值分析上机实习报告学号:姓名:专业:联系电话:任课教师:序 (3)一、必做题 (4)1、问题一 (4)1.1 问题重述 (4)1.2 实验方法介绍 (4)1.3 实验结果 (5)2、问题二 (7)2.1 问题重述 (7)2.2 实验原理 (7)雅各比算法:将系数矩阵A分解为:A=L+U+D,则推到的最后迭代公式为: (8)2.3 实验结果 (8)二、选做题 (10)3、问题三 (10)3.1 问题重述 (10)3.2 实验原理 (10)3.3 实验结果 (11)总结 (11)序伴随着计算机技术的飞速发展,所有的学科都走向定量化和准确化,从而产生了一系列的计算性的学科分支,而数值计算方法就是解决计算问题的桥梁和工具。
数值计算方法,是一种研究并解决数学问题的数值近似解方法,是在计算机上使用的解数学问题的方法。
为了提高计算能力,需要结合计算能力与计算效率,因此,用来解决数值计算的软件因为高效率的计算凸显的十分重要。
数值方法是用来解决数值问题的计算公式,而数值方法的有效性需要根据其方法本身的好坏以及数值本身的好坏来综合判断。
数值计算方法计算的结果大多数都是近似值,但是理论的严密性又要求我们不仅要掌握将基本的算法,还要了解必要的误差分析,以验证计算结果的可靠性。
数值计算一般涉及的计算对象是微积分,线性代数,常微分方程中的数学问题,从而对应解决实际中的工程技术问题。
在借助MA TLAB、JA V A、C++ 和VB软件解决数学模型求解过程中,可以极大的提高计算效率。
本实验采用的是MATLAB软件来解决数值计算问题。
MA TLAB是一种用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境,其对解决矩阵运算、绘制函数/数据图像等有非常高的效率。
本文采用MATLAB对多项式拟合、雅雅格比法与高斯-赛德尔迭代法求解方程组迭代求解,对Runge-Kutta 4阶算法进行编程,并通过实例求解验证了其可行性,使用不同方法对计算进行比较,得出不同方法的收敛性与迭代次数的多少,比较各种方法的精确度和解的收敛速度。
数值代数上机实验报告
数值代数上机实验报告试验项目名称:平方根法与改进平方根法实验内容:先用你熟悉的计算机语言将平方根法和改进平方根法编写成通用的子程序,然后用你编写的程序求解对称正定方程组Ax=b,其中,A=[101 10 1…1 10 11 10]100*100b随机生成,比较计算结果,评论方法优劣。
实验要求:平方根法与改进的平方根的解法步骤;存储单元,变量名称说明;系数矩阵与右端项的生成;结果分析。
实验报告姓名:罗胜利班级:信息与计算科学0802 学号:u200810087 实验一、平方根法与改进平方根法先用你所熟悉的计算机语言将平方根法和改进的平方根法编成通用的子程序,然后用你编写的程序求解对称正定方程组AX=b,其中系数矩阵为40阶Hilbert矩阵,即系数矩阵A的第i行第j列元素为=,向量b的第i个分量为=.平方根法函数程序如下:function [x,b]=pingfanggenfa(A,b)n=size(A);n=n(1);x=A^-1*b; %矩阵求解disp('Matlab自带解即为x');for k=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endend %Cholesky分解for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n); %前代法A=A';for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1); %回代法disp('平方根法的解即为b');endfunction [x]=ave(A,b,n) %用改进平方根法求解Ax=b L=zeros(n,n); %L为n*n矩阵D=diag(n,0); %D为n*n的主对角矩阵S=L*D;for i=1:n %L的主对角元素均为1L(i,i)=1;for i=1:nfor j=1:n %验证A是否为对称正定矩阵if (eig(A)<=0)|(A(i,j)~=A(j,i)) %A的特征值小于0或A非对称时,输出wrong disp('wrong');break;endendendD(1,1)=A(1,1); %将A分解使得A=LDL Tfor i=2:nfor j=1:i-1S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);endD(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');endy=zeros(n,1); % x,y为n*1阶矩阵x=zeros(n,1);for i=1:ny(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i); %通过LDy=b解得y的值endfor i=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n)); %通过L T x=y解得x的值改进平方根法函数程序如下:function b=gaijinpinfanggenfa(A,b)n=size(A);n=n(1);v=zeros(n,1);for j=1:nfor i=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);end %LDL'分解B=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;EndA=tril(A); %得到L和Dfor j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n); %前代法A=D*(A');for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1); %回代法disp('改进平方根法解得的解即为b');end调用函数解题:clear;clc;n=input('请输入矩阵维数:');b=zeros(n,1);A=zeros(n);for i=1:nfor j=1:nA(i,j)=1/(i+j-1);b(i)=b(i)+1/(i+j-1);endend %生成hilbert矩阵[x,b]=pingfanggenfa(A,b) b=gaijinpinfanggenfa(A,b)运行结果:请输入矩阵维数:40Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.570692e-020. > In pingfanggenfa at 4In qiujie at 10Matlab自带解即为x平方根法的解即为bx =1.60358.96850.85621.01950.9375-50.2500-3.0000-16.000024.0000-49.5000-30.000039.000022.0000-64.0000-12.00002.000010.2500-10.5000-1.0000-10.875083.000046.0000-98.0000-69.000068.000021.0000-50.7188-8.7500-8.0000 112.0000 6.0000 -68.7500 22.000044.0000 -28.0000 8.0000 -44.000012.0000b =1.0e+007 *0.0000-0.00000.0001-0.0004-0.00140.0424-0.29801.1419-2.73354.2539-4.30182.7733-1.19890.5406-0.36880.32850.4621-0.25130.05650.0000-0.00510.0071-0.0027-0.0031-0.00190.00090.0002-0.0002-0.00060.00040.0001-0.00020.00010.0000-0.00000.0000-0.0000-0.0000改进平方根法解得的解即为bb =1.0e+024 *0.0000-0.00000.0001-0.0012-0.0954 0.4208 -1.2101 2.0624 -1.0394 -3.3343 6.2567 -0.2463 -7.45942.80303.6990 0.7277 -1.7484 -0.4854 -3.6010 0.2532 5.1862 1.4410 0.8738 -4.5654 1.0422 4.0920 -2.7764 -2.2148 -0.8953 0.3665 4.8967 1.0416 0.1281-1.1902-2.83348.4610-3.6008实验二、利用QR分解解线性方程组:利用QR分解解线性方程组Ax=b,其中A=[16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12];b=[32 26 38 30];求解程序如下:定义house函数:function [v,B]=house(x)n=length(x);y=norm(x,inf);x=x/y;Q=x(2:n)'*x(2:n);v(1)=1;v(2:n)=x(2:n);if n==1B=0;elsea=sqrt(x(1)^2+Q);if x(1)<=0v(1)=x(1)-a;elsev(1)=-Q/(x(1)+a);endB=2*v(1)^2/(Q+v(1)^2);endend进行QR分解:clear;clc;A=[16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12]; b=[32 26 38 30];b=b';x=size(A);m=x(1);n=x(2);d=zeros(n,1);for j=1:n[v,B]=house(A(j:m,j));A(j:m,j:n)=(eye(m-j+1)-B*(v')*v)*A(j:m,j:n); d(j)=B;if j<m< p="">A(j+1:m,j)=v(2:m-j+1);endend %QR分解R=triu(A); %得到R D=A;I=eye(m,n);Q=I;for i=1:nD(i,i)=1;endH=tril(D);M=H';for i=1:nN=I-d(i)*H(1:m,i)*M(i,1:m);Q=Q*N;end %得到Qb=(Q')*b; %Q是正交阵for j=n:-1:2b(j)=b(j)/R(j,j);b(1:j-1)=b(1:j-1)-b(j)*R(1:j-1,j);endb(1)=b(1)/R(1,1); %回带法运行结果如下:R =18.7617 9.8072 15.7769 11.08640 9.9909 9.3358 7.53410 0 5.9945 9.80130 0 0 -0.5126Q =0.8528 -0.4368 -0.2297 -0.17090.2132 0.7916 -0.4594 -0.34170.4264 0.3822 0.2844 0.76890.2132 0.1911 0.8095 -0.5126b=1.000000000000001.000000000000010.9999999999999881.00000000000001实验三、Newton下山法解非线性方程组:3x-cos(yz)-=0,-81+sinz+1.06=0,exp(-xy)+20z+=0;要求满足数值解=满足或.定义所求方程组的函数:Newtonfun.mfunction F = Newtonfun(X)F(1,1)=3*X(1)-cos(X(2)*X(3))-1/2;F(2,1)=X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;F(3,1)=exp(-X(1)*X(2))+20*X(3)+(10*pi-3)/3;End向量求导:Xiangliangqiudao.mfunction J=xiangliangqiudao()syms x y zX=[x,y,z];F=[3*X(1)-cos(X(2)*X(3))-1/2;X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;exp(-X(1)*X(2))+20*X(3)+(10*pi-3)/3];J=jacobian(F,[x y z]);End代值函数:Jacobi.mfunction F=Jacobi(x)F=[ 3,x(3)*sin(x(2)*x(3)), x(2)*sin(x(2)*x(3));2*x(1), -162*x(2)-81/5,cos(x(3));-x(2)/exp(x(1)*x(2)),-x(1)/exp(x(1)*x(2)),20];End方程组求解:format long; %数据表示为双精度型X1=[0,0,0]';eps=10^(-8);k=1;i=1;X2=X1-Jacobi(X1)^(-1)*Newtonfun(X1);while (norm(subs(X2-X1,pi,3.1415926),2)>=eps)&&(norm(Newtonfun(X1),2)>=eps) if norm(Newtonfun(X2),2)<="" p="">X1=X2;B=inv(Jacobi(X2));C=Newtonfun(X2);X2=X2-B*C;i=i+1;elsev=1/(2^k); %引入下山因子X1=X2;B=inv(Jacobi(X2));C=Newtonfun(X2);X2=X2-v*B*C;k=k+1;endendj=i+k-1 %迭代次数X=X2 %输出结果运行结果如下:j =5X =0.500000000000000 -0.000000000000000 -0.523598775598299</m<>。
数值分析上机实验报告
实验报告一题目: (绪论) 非线性方程求解及误差估计摘要:非线性方程的解析解通常很难给出,因此线性方程的数值解法就尤为重要。
本实验采用两种常见的求解方法二分法、Newton 法和改进的Newton 法。
可以节省计算机的计算时间,还能减小不必要的误差。
前言:(目的和意义)掌握二分法与Newton 法的基本原理、应用以及熟练掌握用MATLAB 求函数积分 数学原理:(1)函数的调用格式:quadl(filename,a,b,tol,trace)其中filename 是调用函数名,a 和b 分别为定积分的下限和上限。
用来控制积分精度。
(2)秦九韶算法: S n =a nS k =xS k+1+a k (k=n-1,n-2,...,0), P n (x)=S 0 程序设计:例1.1 计算积分de x xx110利用MATLAB ,下面给出主程序>>g=inline('x.^10.*exp(x-1)'); %定义一个语句函数g(x)=exp(x^10*exp(x-1)) I=quadl(g,0,1) I =0.0098例1.9 秦九韶算法a 0=3,a k=2a k-1+3,Pn(x)=a n x^n+a n-1x^(n-1)+...+a1x+a0求I1=P100(0.5),I2=P150(13)>>x=input('x=');n=input('n=');a=3;for i=1:na=2*a+3;ends=z;b=(a-3)/2;for m=1:100s=x*s+b;b=(b-3)/2;enddisp(s);>>x=0.5n=100600.0000>>x=3n=1004.7039e+078结果分析和讨论:结论:对于二分法,只要能够保证在给定的区间内有根,使能够收敛的,当时收敛的速度和给定的区间有关,二且总体上来说速度比较慢。
暑期实习报告--数值实验
2011年暑期实习报告数值实验姓名: GZ 班级:信息与计算科学09-3班学号:一、 实验目的(1) 《数值计算方法》算法的上机实践。
2011年的春季学习,我们学习了《数值计算方法》这门很重要的课程。
要知道课本上的东西不经过消化是不能成为自己的知识的,不讲理论联系实际,那么我们学到的东西是没有任何意义的。
只有真正去试验,去检验才能成为明白其中的道理。
只有经过怀疑、验证、比较才能明白,什么是好用的,什么是实用的。
为了能够使我们更好的掌握和了解数值计算的算法,为了使我们能将学过的知识真正用于将来的编程设计中,学校和学院老师为我们提供了暑期实践这样一个很好的机会,让大家能亲身感受到算法的奥妙,在短时间内有很大的提高。
(2)Mablab 的学习和使用。
我们在认真领悟算法的同时,通过程序的实现过程,很好的了解并且熟练使用Matlab 也是我们实习的一个重要目的。
一种实用的工具如果不经过练习是不可能被熟练的掌握和使用的,我们如果不经过编程实验,就不能体会到Matlab 强的功能和数学计算上的优势。
Matlab 是三大数学软件之一,有强大的数值计算功能,同时有许多实用的工具箱。
要掌握这些必须要我们多加练习,这次实验正给我们提供了一次了解使用Matlab 编程的机会,在实验中体会Matlab 带给我们的方便。
熟练的掌握Matlab 的用法和工具箱是我们信息专业的同学必须具备的条件。
(3)体会实验的意义和重要性。
提高自己在实验中发现问题、解决问题的能力。
实验的过程是是最重要的,在过程中我们不断的提高自己,发展自己,认识不足,弥补不足才会有真正的提高。
我参加了这次实践,获益良多。
二、 实验内容及结果分析实验一:迭代法求方程根的收敛性内容:为了体会在非线性方程求根的迭代法中,迭代函数和初始迭代值的选取对迭代收敛性的影响。
考虑一个简单的代数方程324100x x +-=,针对上述方程,可以构造多种迭代法,如下列几种迭代格式:①321410n n n x x x x +=--++;②13211(10)2n n x x +=-;③12110()4n x x +=+;④321241038n n n n n nx x x x x x ++-=-+。
数值分析2024上机实验报告
数值分析2024上机实验报告数值分析是计算数学的一个重要分支,它研究如何用数值方法来解决数学问题。
在数值分析的学习过程中,学生需要通过上机实验来巩固理论知识,并学会使用相应的数值方法来解决实际问题。
本篇报告将详细介绍2024年度数值分析上机实验的内容和结果。
一、实验内容2024年度数值分析上机实验分为四个部分,分别是:方程求根、插值与拟合、数值积分和常微分方程的数值解。
1.方程求根这部分实验要求使用数值方法求解给定的非线性方程的根。
常见的数值方法有二分法、牛顿法、割线法等。
在实验过程中,我们需要熟悉这些数值方法的原理和实现步骤,并对不同方法的收敛性进行分析和比较。
2.插值与拟合这部分实验要求使用插值和拟合方法对给定的一组数据进行拟合。
插值方法包括拉格朗日插值、牛顿插值等;拟合方法包括最小二乘拟合、多项式拟合等。
在实验中,我们需要熟悉插值和拟合方法的原理和实现步骤,并对不同方法的精度和稳定性进行比较。
3.数值积分这部分实验要求使用数值方法计算给定函数的积分。
常见的数值积分方法有梯形法则、辛普森法则、龙贝格积分等。
在实验过程中,我们需要熟悉这些数值积分方法的原理和实现步骤,并对不同方法的精度和效率进行比较。
4.常微分方程的数值解这部分实验要求使用数值方法求解给定的常微分方程初值问题。
常见的数值方法有欧拉法、改进的欧拉法、四阶龙格-库塔法等。
在实验中,我们需要熟悉这些数值解方法的原理和实现步骤,并对不同方法的精度和稳定性进行比较。
二、实验结果在完成2024年度数值分析上机实验后,我们得到了以下实验结果:1.方程求根我们实现了二分法、牛顿法和割线法,并对比了它们的收敛速度和稳定性。
结果表明,割线法的收敛速度最快,但在一些情况下可能会出现振荡;二分法和牛顿法的收敛速度相对较慢,但稳定性较好。
2.插值与拟合我们实现了拉格朗日插值和最小二乘拟合,并对比了它们的拟合效果和精度。
结果表明,拉格朗日插值在小区间上拟合效果较好,但在大区间上可能出现振荡;最小二乘拟合在整体上拟合效果较好,但可能出现过拟合。
数值分析上机实习报告
指导教师:姓名:学号:专业:联系电话:上海交通大学目录序言 (3)实验课题(一) 雅可比迭代法和高斯-塞得尔迭代法的收敛性和收敛速度 (4)数值分析 (6)实验课题(二) 松弛因子对SOR法收敛速度的影响 (6)数值分析 (12)总结 (13)附录(程序清单) (14)1.雅可比迭代法和高斯-塞得尔迭代法的收敛性和收敛速度 (14)雅可比迭代法: (14)高斯-塞得尔迭代法: (16)2.松弛因子对SOR法收敛速度的影响 (18)松弛法(SOR) (18)序言随着科学技术的发展,提出了大量复杂的数值计算问题,在实际解决这些计算问题的长期过程中,形成了计算方法这门学科,专门研究各种数学问题的数值解法(近似解法),包括方法的构造和求解过程的误差分析,是一门内容丰富,有自身理论体系的实用性很强的学科。
解决工程问题,往往需要处理很多数学模型,这就要花费大量的人力和时间,但是还有不少数学模型无法用解析法得到解。
使用数值方法并利用计算机,就可以克服这些困难。
事实上,科学计算已经与理论分析、科学实验成为平行的研究和解决科技问题的科学手段,经常被科技工作者所采用。
作为科学计算的核心内容——数值分析(数值计算方法),已逐渐成为广大科技工作者必备的基本知识并越来越被人重视。
由于数值方法是解数值问题的系列计算公式,所以数值方法是否有效,不但与方法本身的好坏有关,而且与数值问题本身的好坏也有关,因此,研究数值方法时,不但需要研究数值方法的好坏,即数值稳定性问题,而且还需要研究数值问题本身的好坏,即数值问题的性态,以及它们的判别问题。
数值计算的绝大部分方法都具有近似性,而其理论又具有严密的科学性,方法的近似值正是建立在理论的严密性基础上,根据计算方法的这一特点。
因此不仅要求掌握和使用算法,还要重视必要的误差分析,以保证计算结果的可靠性。
数值计算还具有应用性强的特点,计算方法的绝大部分方法如求微分方程近似解,求积分近似值,求解超越方程,解线性方程组等都具有较强的实用性,而插值法,最小二乘法,样条函数等也都是工程技术领域中常用的,有实际应用价值的方法。
(完整word版)计算方法A上机实验报告
计算方法A上机实验报告姓名:苏福班级:硕4020 学号:3114161019一、上机练习目的1)复习和巩固数值计算方法的基本数学模型,全面掌握运用计算机进行数值计算的具体过程及相关问题。
2)利用计算机语言独立编写、调试数值计算方法程序,培养学生利用计算机和所学理论知识分析解决实际问题的能力。
二、上机练习任务1)利用计算机语言编写并调试一系列数值方法计算通用程序,并能正确计算给定题目,掌握调试技能。
2)掌握文件使用编程技能,如文件的各类操作,数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。
3)写出上机练习报告。
三、上机题目1. 共轭梯度法求解线性方程组。
(第三章)2. 三次样条插值(第四章)3. 龙贝格积分(第六章)4. 四阶龙格-库塔法求解常微分方程的初值问题四、上机报告题目1:共轭梯度法求解线性方程组1.算法原理共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小值的问题。
从任意给定的初始点出发,沿一组关于矩阵A共轭的方向进行线性搜索,在无舍入误差的假定下,最多迭代n 次(其中n 为矩阵A 的阶数),就可求得二次函数的极小值,也就求得了线性方程组Ax b =的解。
定理:设A 是n 阶对称正定矩阵,则x *是方程组Ax b =的解得充分必要条件是x *是二次函数1()2TT f x x Ax b x =-的极小点,即 ()()min nx R Ax b f x f x **∈=⇔=共轭梯度法的计算公式:(0)(0)(0)()()()()(1)()()(1)(1)(1)()()()(1)(1)()k T k k k T k k k k k k k k T k k k T k k k k k d r b Ax r d d Ad xx d r b Ax r Ad d Ad d r d ααββ++++++⎧==-⎪⎪=⎪⎪=+⎪⎨=-⎪⎪⎪=-⎪⎪=+⎩2. 程序框图(1)编写共轭梯度法求解对称正定矩阵的线性方程组见附录(myge.m):function x=myge(A,b)输入对称正定矩阵及对应的列向量,初始向量设为0,精度取为810 。
计算方法上机实验
1.拉格朗日插值多项式,用于离散数据的拟合#include <stdio.h>#include <conio.h>#include <alloc.h>float lagrange(float *x,float *y,float xx,int n) /*拉格朗日插值算法*/{ int i,j;float *a,yy=0.0; /*a作为临时变量,记录拉格朗日插值多项式*/a=(float *)malloc(n*sizeof(float));for(i=0;i<=n-1;i++){ a[i]=y[i];for(j=0;j<=n-1;j++)if(j!=i) a[i]*=(xx-x[j])/(x[i]-x[j]);yy+=a[i];}free(a);return yy;}main(){ int i,n;float x[20],y[20],xx,yy;printf("Input n:");scanf("%d",&n);if(n>=20) {printf("Error!The value of n must in (0,20)."); getch();return 1;} if(n<=0) {printf("Error! The value of n must in (0,20)."); getch(); return 1;} for(i=0;i<=n-1;i++){ printf("x[%d]:",i);scanf("%f",&x[i]);}printf("\n");for(i=0;i<=n-1;i++){ printf("y[%d]:",i);scanf("%f",&y[i]);}printf("\n");printf("Input xx:");scanf("%f",&xx);yy=lagrange(x,y,xx,n);printf("x=%f,y=%f\n",xx,yy);getch();}2.牛顿插值多项式,用于离散数据的拟合#include <stdio.h>#include <conio.h>#include <alloc.h>void difference(float *x,float *y,int n){ float *f;int k,i;f=(float *)malloc(n*sizeof(float));for(k=1;k<=n;k++){ f[0]=y[k];for(i=0;i<k;i++)f[i+1]=(f[i]-y[i])/(x[k]-x[i]);y[k]=f[k];}return;}main(){ int i,n;float x[20],y[20],xx,yy;printf("Input n:");scanf("%d",&n);if(n>=20) {printf("Error! The value of n must in (0,20)."); getch(); return 1;} if(n<=0) {printf("Error! The value of n must in (0,20).");getch(); return 1;} for(i=0;i<=n-1;i++){ printf("x[%d]:",i);scanf("%f",&x[i]);}printf("\n");for(i=0;i<=n-1;i++){ printf("y[%d]:",i);scanf("%f",&y[i]);}printf("\n");difference(x,(float *)y,n);printf("Input xx:");scanf("%f",&xx);yy=y[20];for(i=n-1;i>=0;i--) yy=yy*(xx-x[i])+y[i];printf("NewtonInter(%f)=%f",xx,yy);getch();}3.高斯列主元消去法,求解其次线性方程组第一种#include<stdio.h>#include <math.h>#define N 20int main(){ int n,i,j,k;int mi,tmp,mx;float a[N][N],b[N],x[N];printf("\nInput n:");scanf("%d",&n);if(n>N){ printf("The input n should in(0,N)!\n");getch();return 1;}if(n<=0){ printf("The input n should in(0,N)!\n");getch();return 1;}printf("Now input a(i,j),i,j=0...%d:\n",n-1); for(i=0;i<n;i++){ for(j=0;j<n;j++)scanf("%f",&a[i][j]);}printf("Now input b(i),i,j=0...%d:\n",n-1);for(i=0;i<n;i++)scanf("%f",&b[i]);for(i=0;i<n-2;i++){ for(j=i+1,mi=i,mx=fabs(a[i][j]);j<n-1;j++) if(fabs(a[j][i])>mx){ mi=j;mx=fabs(a[j][i]);}if(i<mi){ tmp=b[i];b[i]=b[mi];b[mi]=tmp;for(j=i;j<n;j++){ tmp=a[i][j];a[i][j]=a[mi][j];a[mi][j]=tmp;}}for(j=i+1;j<n;j++){ tmp=-a[j][i]/a[i][i];b[j]+=b[i]*tmp;for(k=i;k<n;k++)a[j][k]+=a[i][k]*tmp;}}x[n-1]=b[n-1]/a[n-1][n-1];for(i=n-2;i>=0;i--){ x[i]=b[i];for(j=i+1;j<n;j++)x[i]-=a[i][j]*x[j];x[i]/=a[i][i];}for(i=0;i<n;i++)printf("Answer:\n x[%d]=%f\n",i,x[i]); getch();return 0;}第二种#include<math.h>#include<stdio.h>#define NUMBER 20#define Esc 0x1b#define Enter 0x0dfloat A[NUMBER][NUMBER+1] ,ark;int flag,n;exchange(int r,int k);float max(int k);message();main(){float x[NUMBER];int r,k,i,j;char celect;clrscr();printf("\n\nUse Gauss.");printf("\n\n1.Jie please press Enter."); printf("\n\n2.Exit press Esc.");celect=getch();if(celect==Esc)exit(0);printf("\n\n input n=");scanf("%d",&n);printf(" \n\nInput matrix A and B:"); for(i=1;i<=n;i++){printf("\n\nInput a%d1--a%d%d and b%d:",i,i,n,i);for(j=1;j<=n+1;j++) scanf("%f",&A[i][j]);}for(k=1;k<=n-1;k++){ark=max(k);if(ark==0){printf("\n\nIt's wrong!");message();}else if(flag!=k)exchange(flag,k);for(i=k+1;i<=n;i++)for(j=k+1;j<=n+1;j++)A[i][j]=A[i][j]-A[k][j]*A[i][k]/A[k][k];}x[n]=A[n][n+1]/A[n][n];for( k=n-1;k>=1;k--){float me=0;for(j=k+1;j<=n;j++){me=me+A[k][j]*x[j];}x[k]=(A[k][n+1]-me)/A[k][k];}for(i=1;i<=n;i++){printf(" \n\nx%d=%f",i,x[i]);}message();}exchange(int r,int k){int i;for(i=1;i<=n+1;i++)A[0][i]=A[r][i];for(i=1;i<=n+1;i++)A[r][i]=A[k][i];for(i=1;i<=n+1;i++)A[k][i]=A[0][i];}float max(int k){int i;float temp=0;for(i=k;i<=n;i++)if(fabs(A[i][k])>temp){temp=fabs(A[i][k]);flag=i;}return temp;}message(){printf("\n\n Go on Enter ,Exit press Esc!");switch(getch()){case Enter: main();case Esc: exit(0);default:{printf("\n\nInput error!");message();} }}4.龙贝格求积公式,求解定积分#include<stdio.h>#include<math.h>#define f(x) (sin(x)/x)#define N 20#define MAX 20#define a 2#define b 4#define e 0.00001float LBG(float p,float q,int n){ int i;float sum=0,h=(q-p)/n;for (i=1;i<n;i++)sum+=f(p+i*h);sum+=(f(p)+f(q))/2;return(h*sum);}void main(){ int i;int n=N,m=0;float T[MAX+1][2];T[0][1]=LBG(a,b,n);n*=2;for(m=1;m<MAX;m++){ for(i=0;i<m;i++)T[i][0]=T[i][1];T[0][1]=LBG(a,b,n);n*=2;for(i=1;i<=m;i++)T[i][1]=T[i-1][1]+(T[i-1][1]-T[i-1][0])/(pow(2,2*m)-1);if((T[m-1][1]<T[m][1]+e)&&(T[m-1][1]>T[m][1]-e)){ printf("Answer=%f\n",T[m][1]); getch();return ;}}}5.牛顿迭代公式,求方程的近似解#include<stdio.h>#include<math.h>#include<conio.h>#define N 100#define PS 1e-5#define TA 1e-5float Newton(float (*f)(float),float(*f1)(float),float x0 ) { float x1,d=0;int k=0;do{ x1= x0-f(x0)/f1(x0);if((k++>N)||(fabs(f1(x1))<PS)){ printf("\nFailed!");getch();exit();}d=(fabs(x1)<1?x1-x0:(x1-x0)/x1);x0=x1;printf("x(%d)=%f\n",k,x0);}while((fabs(d))>PS&&fabs(f(x1))>TA) ;return x1;}float f(float x){ return x*x*x+x*x-3*x-3; }float f1(float x){ return 3.0*x*x+2*x-3; }void main(){ float f(float);float f1(float);float x0,y0;printf("Input x0: ");scanf("%f",&x0);printf("x(0)=%f\n",x0);y0=Newton(f,f1,x0);printf("\nThe root is x=%f\n",y0); getch();}。
数值计算方法上机实验报告
数值计算方法上机实验报告
一、实验目的
本次实验的主要目的是熟悉和掌握数值计算方法,学习梯度下降法的
原理和实际应用,熟悉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库的使用方法。
数值分析上机实验报告
数值分析上机实验理学院11级统计01班41108030125鲁庆实验报告一一.实验名称误差与误差估计二.实验目的掌握数值运算的误差估计方法三.数学原理 1.绝对误差(*)e x设某一量的准确值为x ,近似值为x*,则x*与x 之差叫做近似值x*的绝对误差(简称误差),记为*(*)*e e x x x ==- 2.绝对误差限适当小的正数,使|(*)||*|*e x x x ε=-≤则称*ε为近似值 x * 的绝对误差限。
(有时用*x x ε*=±表示近似值x *的精度或准确值的所在范围。
3.相对误差(*)r e x绝对误差与准确值之比*(*)*(*),0r r e x x xe e x x x x-===≠称为x *的相对 误差。
4.相对误差限(*)r x ε若指定一个适当小的正数 (*)r x ε,使|(*)||(*)|(*)||r r e x e x x x ε=≤则称(*)r x ε为近似值 x *的相对误差限。
5.有效数字若近似值x*的绝对误差限是某一位的半个单位,该位到x*的第一位非零数字一共有n 位,则称近似值x*有n 位有效数字,或说x*精确到该位。
6.绝对误差的运算:)()()(2121x x x x εεε+=± )()()(122121x x x x x x εεε+≈22122121+=x x x x x x x )()()(εεε (f(x))()(x)f x εε'≈四.实验内容1. 计算I n=e 1-⎰10nxe x 2dx (n=0,1,...)并估计误差。
解: >> I0 = exp(-1)*quad('(x.^0).*exp(x.^2)',0,1,10^(-10));>> vpa(I0,10) ans =.5380795069>> I1= exp(-1)*quad('(x.^1).*exp(x.^2)',0,1,10^(-10)); >> vpa(I1,10) ans =.3160602794>> I2 = exp(-1)*quad('(x.^2).*exp(x.^2)',0,1,10^(-10)); >> vpa(I2,10) ans =.2309602465>> I3 = exp(-1)*quad('(x.^3).*exp(x.^2)',0,1,10^(-10)); >> vpa(I3,10) ans =.1839397206>> I4 = exp(-1)*quad('(x.^4).*exp(x.^2)',0,1,10^(-10)); >> vpa(I4,10) ans =.1535596302>> I5 = exp(-1)*quad('(x.^5).*exp(x.^2)',0,1,10^(-10)); >> vpa(I5,10) ans =.1321205588>> I6 = exp(-1)*quad('(x.^6).*exp(x.^2)',0,1,10^(-10)); >> vpa(I6,10) ans =.1161009245>> I7 = exp(-1)*quad('(x.^7).*exp(x.^2)',0,1,10^(-10)); >> vpa(I7,10) ans =.1036383235>> I8 = exp(-1)*quad('(x.^8).*exp(x.^2)',0,1,10^(-10)); >> vpa(I8,10) ans =.9364676413e-1>> I9 = exp(-1)*quad('(x.^9).*exp(x.^2)',0,1,10^(-10)); >> vpa(I9,10) ans =.8544670595e-1 2.计算x255的值。
计算方法与实习上机实验报告
计算方法与实习上机实验报告一、引言本文旨在介绍和展示我们在“计算方法”课程中的实习上机实验环节所完成的一些关键任务和所取得的成果。
该实验课程的目标是让我们更深入地理解和应用各种计算方法,并在实际操作中提高我们的编程和问题解决能力。
二、实验内容与目标实验的主要内容是利用各种计算方法解决实际数学问题。
我们被要求使用编程语言(如Python或Java)来实现和解决这些问题。
这些问题包括使用牛顿法求解平方根,使用蒙特卡洛方法计算圆周率,以及使用最优化方法求解函数的最小值等。
实验的目标不仅是让我们掌握计算方法的基本理论,更是要让我们能够在实际操作中运用这些方法。
我们需要在实习过程中,通过与同伴们合作,共同解决问题,提高我们的团队合作能力和问题解决能力。
三、实验过程与问题解决策略在实验过程中,我们遇到了许多问题,如编程错误、理解困难和时间压力等。
我们通过相互讨论、查阅资料和寻求教师帮助等方式,成功地解决了这些问题。
例如,在实现牛顿法求解平方根时,我们一开始对导数的计算和理解出现了一些错误。
但我们通过查阅相关资料和讨论,最终理解了导数的正确计算方法,并成功地实现了牛顿法。
四、实验结果与结论通过这次实习上机实验,我们不仅深入理解了计算方法的基本理论,还在实际操作中提高了我们的编程和问题解决能力。
我们的成果包括编写出了能有效求解平方根、计算圆周率和求解函数最小值的程序。
这次实习上机实验非常成功。
我们的团队不仅在理论学习和实践操作上取得了显著的进步,还在团队合作和问题解决方面积累了宝贵的经验。
这次实验使我们对计算方法有了更深的理解和认识,也提高了我们的编程技能和解决问题的能力。
五、反思与展望回顾这次实验,我们意识到在实验过程中,我们需要更好地管理我们的时间和压力。
在解决问题时,我们需要更有效地利用我们的知识和资源。
在未来,我们希望能够更加熟练地运用计算方法,并能够更有效地解决问题。
我们也希望能够将所学的计算方法应用到更广泛的领域中,如数据分析、科学研究和工业生产等。
计算方法与实习上机报告
计算方法与实习——上机报告学院:电子工程学院学号:姓名:刘波2015.1.4计算方法与实习上机报告习题一:1 舍入误差及稳定性一、实验目的(1)通过上机编程,复习巩固以前所学程序设计语言及上机操作指令;(2)通过上机计算,了解舍入误差所引起的数值不稳定性二、实验内容1、用两种不同的顺序计算1000021n n -=∑,分析其误差的变化 2、已知连分数()101223//(.../)n n a f b b a b a a b =++++,利用下面的算法计算f : 11,i n n i i i a d b d b d ++==+(1,2,...,0)i n n =-- 0f d = 写一程序,读入011,,,...,,,...,,n n n b b b a a 计算并打印f 3、给出一个有效的算法和一个无效的算法计算积分1041nn x y dx x =+⎰ (0,1,...,10)n = 4、设2211N N j S j ==-∑,已知其精确值为1311221N N ⎛⎫-- ⎪+⎝⎭ (1)编制按从大到小的顺序计算N S 的程序 (2)编制按从小到大的顺序计算N S 的程序(3)按两种顺序分别计算10001000030000,,,S S S 并指出有效位数三、实验步骤、程序设计、实验结果及分析1、用两种不同的顺序计算1000021n n -=∑,分析其误差的变化 (1)实验步骤:分别从1~10000和从10000~1两种顺序进行计算,应包含的头文件有stdio.h 和math.h(2)程序设计:a.顺序计算#include<stdio.h>#include<math.h>void main(){double sum=0;int n=1;while(1){sum=sum+(1/pow(n,2));if(n%1000==0)printf("sun[%d]=%-30f",n,sum);if(n>=10000)break;n++;}printf("sum[%d]=%f\n",n,sum);}b.逆序计算#include<stdio.h>#include<math.h>void main(){double sum=0;int n=10000;while(n!=0){sum=sum+(1/pow(n,2));if(n%200==0)printf("sum[%d]=%-10f",n,sum);if(n<1)break;n--;}printf("sum[%d]=%f\n",n,sum);}(3)实验结果及分析:程序运行结果:a.顺序计算b.逆序计算结果分析:两种不同顺序计算结果是一样的,顺序计算误差从一开始就很小,而逆序计算误差最开始十分大,后来结果正确。
计算方法上机实验报告
. 《计算方法》上机实验报告班级:XXXXXX小组成员:XXXXXXXXXXXXXXXXXXXXXXXXXXXX任课教师:XXX二〇一八年五月二十五日前言通过进行多次的上机实验,我们结合课本上的内容以及老师对我们的指导,能够较为熟练地掌握Newton 迭代法、Jacobi 迭代法、Gauss-Seidel 迭代法、Newton 插值法、Lagrange 插值法和Gauss 求积公式等六种算法的原理和使用方法,并参考课本例题进行了MATLAB 程序的编写。
以下为本次上机实验报告,按照实验内容共分为六部分。
实验一:一、实验名称及题目: Newton 迭代法例2.7(P38):应用Newton 迭代法求在附近的数值解,并使其满足.二、解题思路:设'x 是0)(=x f 的根,选取0x 作为'x 初始近似值,过点())(,00x f x 做曲线)(x f y =的切线L ,L 的方程为))((')(000x x x f x f y -+=,求出L 与x 轴交点的横坐标)(')(0001x f x f x x -=,称1x 为'x 的一次近似值,过点))(,(11x f x 做曲线)(x f y =的切线,求该切线与x 轴的横坐标)(')(1112x f x f x x -=称2x 为'x 的二次近似值,重复以上过程,得'x 的近似值序列{}n x ,把)(')(1n n n n x f x f x x -=+称为'x 的1+n 次近似值,这种求解方法就是牛顿迭代法。
三、Matlab 程序代码:function newton_iteration(x0,tol) syms z %定义自变量 format long %定义精度 f=z*z*z-z-1; f1=diff(f);%求导 y=subs(f,z,x0);y1=subs(f1,z,x0);%向函数中代值 x1=x0-y/y1; k=1;while abs(x1-x0)>=tol x0=x1;y=subs(f,z,x0); y1=subs(f1,z,x0); x1=x0-y/y1;k=k+1; endx=double(x1) K四、运行结果:实验二:一、实验名称及题目:Jacobi 迭代法例3.7(P74):试利用Jacobi 迭代公式求解方程组要求数值解为方程组的精确解. 二、解题思路:首先将方程组中的系数矩阵A 分解成三部分,即:U D L A ++=,D 为对角阵,L 为下三角矩阵,U 为上三角矩阵。
实用数值计算方法上机实验报告
实用数值计算方法上机实验报告学院:化学工程学院姓名:**专业:工业催化学号: **********1. 问题来源某公司饲养实验用的动物以供出售,已知这些动物的生长对饲料中3种营养成分(蛋白质,矿物质和维生素)特别敏感,每个动物每周至少需要蛋白质60g,矿物质3g,维生素8mg,该公司能买到5种不同的饲料,每种饲料1kg 所含各种营养成分和成本如表1所示,如果每个小动物每周食用饲料不超过52kg,求既满足动物生长需要,又能使总成本最低的饲料配方。
数学模型 设需要饲料A1,A2,A3,A4,A5分别为x1,x2,x3,x4,x5(单位kg )12345min 0.20.70.40.30.5S x x x x x =++++1234512345123451234512350.3x +2x +x +0.6x +1.8x 600.1x +0.05x +0.02x +0.2x +0.05x 3.0.05x +0.1x +0.02x +0.2x +0.08x 8x +x +x +x +x 52,,,4,0s t x x x x x ≥⎧⎪≥⎪⎪≥⎨⎪≤⎪⎪≥⎩在LINGO 的MODEL 窗口内输入如下模型:Min =0.2*x1+0.7*x2+0.4*x3+0.3*x4+0.5*x5; 0.3*x1+2*x2+x3+0.6*x4+1.8*x5>60;0.1*x1+0.05*x2+0.02*x3+0.2*x4+0.05*x5>3; 0.05*x1+0.1*x2+0.02*x3+0.2*x4+0.08*x5>8; x1+x2+x3+x4+x5<52; end求解输出结果如下:Global optimal solution found.Objective value: 22.40000Infeasibilities: 0.000000Total solver iterations: 3Variable Value Reduced Cost X1 0.000000 0.7000000 X2 12.00000 0.000000 X3 0.000000 0.6166667 X4 30.00000 0.000000 X5 10.00000 0.000000 Row Slack or Surplus Dual Price1 22.40000 -1.0000002 0.000000 -0.58333333 4.100000 0.0000004 0.000000 -4.1666675 0.000000 0.8833333结果分析:因此每周每个动物的配料为饲料A2,A4,A5分别为12kg,30kg,10kg,可使得成本达到最低,最低成本为22.4元。
数值分析上机实践报告
数值分析上机实践报告班级:计算机1002姓名:陈斯琪学号:20102686课题三A . 实验题目:线性方程组的迭代法B . 实验要求(1) 应用迭代法求解线性方程组,并与直接法作比较;(2) 分别对不同精度要求,如5-4-3-10,10,10=ε,利用所需迭代次数体会该迭代法的收敛快慢;(3) 对方程组(2),(3)使用SOR 方法时,选取松弛因子=0.8,0.9,1,1.1,1.2等,试观察对算法收敛性的影响,并找出你所选用松弛因子的最佳值; (4) 编制出各种迭代法的程序并给出计算结果。
C . 目的和意义(1) 通过上机了解迭代法求解线性方程组的特点;掌握求解线性方程组的各类迭代法;(2) 体会上机计算时,终止准则‖X^(k+1)-X^k ‖∞<ε,对控制迭代精度的有效性; (3) 体会初始值和松弛因子的选择,对迭代收敛速度的影响 D . 实验方程组(1)线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡1-421534100368-24-3-81-012029137-2621-234179-11-1003524-31-23-6217758-6233-761-62911-31-512-301-231-2-2010563-5-6000121-3-20416084-0484⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡10987654321x x x x x x x x x x =⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-2119381346323125精确解Tx )2,1,1,3,0,2,1,0,1,1(*--=.(2) 对称正定线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡45152211236601924-3-360024-3-36014110-3-5211144-3-310-4221-8-13-4-1-612-53-8-1141-2312-1-204204-2004204-2487654321x x x x x x x x精确解T*)2,0,1,1,2,0,1,1(--=x .(3)三对角线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡554141262135741-000000001-0041-0000001-41-0000001-41-0000001-41-0000001-41-0000001-41-0000001-41-0000001-400000001-000000001-410987654321x x xx x x x x x x精确解Tx )1,1,0,3,2,1,0,3,1,2(*---=.E . 实验程序代码及截图(1) 应用Jacobi 迭代法求解方程组代码如下: #include<iostream.h> #include<math.h>#define N 10 //十阶矩阵 staticdoubleA[N][N]={4,2,-3,-1,2,1,0,0,0,0,8,6,-5,-3,6,5,0,1,0,0,4,2,-2,-1,3,2,-1,0,3,1,0,-2,1,5,-1,3,-1,1,9,4,-4,2,6,-1,6,7,-3,3,2,3,8,6,-8,5,7,17,2,6,-3,5,0,2,-1,3,-4,2,5,3,0,1,16,10,-11,-9,17,34,2,-1,2,2,4,6,2,-7,13,9,2,0,12,4,0,0,-1,8,-3,-24,-8,6,3,-1};//方程组左侧系数矩阵 static double B[N]={5,12,3,2,3,46,13,38,19,-21}; //右侧值static double Y[N]; //输出比较项static double Y[N];static double X[N]; //输出项static double G[N]; //X = BX' + G的G矩阵int i,j,k; //计数器double eps;int M=100;bool distance(){ //求两输出项的差的范数是否满足精度要求double temp=0;for (i=0;i<N;i++){temp=temp+fabs(X[i]-Y[i]);}if (temp>eps)return false;elsereturn true; //满足精度要求则结束程序}void main(){cout<<"最大迭代次数为100次"<<endl;cout<<"你希望的精度是多少?"<<endl;cout<<"eps=";cin>>eps;//形成迭代矩阵B,存放到A中for (i=0;i<N;i++){if (fabs(A[i][i])<eps){cout <<"打印失败"<<endl;return;}double T=A[i][i];for (j=0;j<N;j++){A[i][j]=-A[i][j]/T;}A[i][i] = 0;G[i]=B[i]/T;}int counter=0;while (counter<M){//迭代for (i=0;i<N;i++){double temp=0;for (j=0;j<N;j++){temp=temp+A[i][j]*Y[j];}X[i]=G[i]+temp;}if (distance()==true)break;else{//交换X,Y向量;for(i=0;i<N;i++){Y[i]=X[i];}}counter++;}//打印Xcout << "迭代次数为:"<<counter<<"次。
数值计算方法上机实验报告
数值计算方法上机实验报告实验目的:复习和巩固数值计算方法的基本数学模型,全面掌握运用计算机进行数值计算的具体过程及相关问题。
利用计算机语言独立编写、调试数值计算方法程序,培养学生利用计算机和所学理论知识分析解决实际问题的能力。
上机练习任务:利用计算机基本C 语言编写并调试一系列数值方法计算通用程序,并能正确计算给定题目,掌握调试技能。
掌握文件使用编程技能,如文件的各类操作,数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。
一、各算法的算法原理及计算机程序框图1. 列主元高斯消去法算法原理:高斯消去法是利用现行方程组初等变换中的一种变换,即用一个不为零的数乘一个方程后加只另一个方程,使方程组变成同解的上三角方程组,然后再自下而上对上三角方程组求解。
列选住院是当高斯消元到第k 步时,从k 列的kk a 以下(包括kk a )的各元素中选出绝对值最大的,然后通过行交换将其交换到kk a 的位置上。
交换系数矩阵中的两行(包括常数项),只相当于两个方程的位置交换了,因此,列选主元不影响求解的结果。
●源程序:#define N 200#include "stdio.h"#include "math.h"FILE *fp1,*fp2;void LZ(){int n,i,j,k=0,l;double d,t,t1;static double x[N],a[N][N];fp1=fopen("a1.txt","r");fp2=fopen("b1.txt","w");fscanf(fp1,"%d",&n);for(i=0;i<n;++i)for(j=0;j<=n;++j){fscanf(fp1,"%lf",&a[i][j]);}{d=a[k][k];l=k;i=k+1;do{if(fabs(a[i][k])>fabs(d)) /*选主元*/{d=a[i][k];l=i;}i++;}while(i<n);if(d==0){printf("\n输入矩阵有误!\n");}else{ /*换行*/if(l!=k){for(j=k;j<=n;j++){t=a[l][j];a[l][j]=a[k][j];a[k][j]=t;}}}for(j=k+1;j<=n;j++) /*正消*/ a[k][j]/=a[k][k];for(i=k+1;i<n;i++)for(j=k+1;j<=n;j++)a[i][j]-=a[i][k]*a[k][j];k++;}while(k<n);if(k!=0){for(i=n-1;i>=0;i--) /*回代*/ {t1=0;for(j=i+1;j<n;j++)t1+=a[i][j]*x[j];x[i]=a[i][n]-t1;}for(i=0;i<n;i++)fprintf(fp2,"\n 方程组的根为x[%d]=%lf",i+1,x[i]); fclose(fp1); fclose(fp2); }main() { LZ(); }● 具体算例及求解结果:用列选主元法求解下列线性方程组⎪⎩⎪⎨⎧=++=++=-+28x x 23x 2232832321321321x x x x x x 输入3 输出结果:方程组的根为x[1]=6.0000001 2 -3 8 方程组的根为x[2]=4.000000 2 1 3 22 方程组的根为x[3]=2.000000 3 2 1 28● 输入变量、输出变量说明:输入变量:ij a 系数矩阵元素,i b 常向量元素 输出变量:12,,n b b b 解向量元素2. 杜里特尔分解法解线性方程● 算法原理:求解线性方程组Ax b =时,当对A 进行杜里特尔分解,则等价于求解LUx b =,这时可归结为利用递推计算相继求解两个三角形(系数矩阵为三角矩阵)方程组,用顺代,由Ly b =求出y ,再利用回带,由Ux y =求出x 。
(完整版)哈工大-数值分析上机实验报告
实验报告一题目:非线性方程求解摘要:非线性方程的解析解通常很难给出,因此线性方程的数值解法就尤为重要。
本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。
前言:(目的和意义)掌握二分法与Newton法的基本原理和应用。
数学原理:对于一个非线性方程的数值解法很多。
在此介绍两种最常见的方法:二分法和Newton法。
对于二分法,其数学实质就是说对于给定的待求解的方程f(x),其在[a,b]上连续,f(a)f(b)<0,且f(x)在[a,b]内仅有一个实根x*,取区间中点c,若,则c恰为其根,否则根据f(a)f(c)<0是否成立判断根在区间[a,c]和[c,b]中的哪一个,从而得出新区间,仍称为[a,b]。
重复运行计算,直至满足精度为止。
这就是二分法的计算思想。
Newton法通常预先要给出一个猜测初值x0,然后根据其迭代公式产生逼近解x*的迭代数列{x k},这就是Newton法的思想。
当x0接近x*时收敛很快,但是当x0选择不好时,可能会发散,因此初值的选取很重要。
另外,若将该迭代公式改进为其中r为要求的方程的根的重数,这就是改进的Newton法,当求解已知重数的方程的根时,在同种条件下其收敛速度要比Newton法快的多。
程序设计:本实验采用Matlab的M文件编写。
其中待求解的方程写成function的方式,如下function y=f(x);y=-x*x-sin(x);写成如上形式即可,下面给出主程序。
二分法源程序:clear%%%给定求解区间b=1.5;a=0;%%%误差R=1;k=0;%迭代次数初值while (R>5e-6) ;c=(a+b)/2;if f12(a)*f12(c)>0;a=c;elseb=c;endR=b-a;%求出误差k=k+1;endx=c%给出解Newton法及改进的Newton法源程序:clear%%%% 输入函数f=input('请输入需要求解函数>>','s')%%%求解f(x)的导数df=diff(f);%%%改进常数或重根数miu=2;%%%初始值x0x0=input('input initial value x0>>');k=0;%迭代次数max=100;%最大迭代次数R=eval(subs(f,'x0','x'));%求解f(x0),以确定初值x0时否就是解while (abs(R)>1e-8)x1=x0-miu*eval(subs(f,'x0','x'))/eval(subs(df,'x0','x'));R=x1-x0;x0=x1;k=k+1;if (eval(subs(f,'x0','x'))<1e-10);breakendif k>max;%如果迭代次数大于给定值,认为迭代不收敛,重新输入初值ss=input('maybe result is error,choose a new x0,y/n?>>','s');if strcmp(ss,'y')x0=input('input initial value x0>>');k=0;elsebreakendendendk;%给出迭代次数x=x0;%给出解结果分析和讨论:1.用二分法计算方程在[1,2]内的根。
数值计算方法上机实验报告
(2)计算机程序框图:(见下页)
(3)输入变量、输出变量说明:
输入变量: 系数矩阵元素, 常向量元素
称为改进欧拉公式。
(2)计算机程序框图:(见下页)
(3)输入变量、输出变量说明:
输入变量: 处置点, 区间长度, 计算次数
输出变量: 初值问题的数值解法结果
(4)具体算例及求解结果:
例:求解初值问题(课本P242例7-2)
求解结果:
0.1
1.095909
1.095909
0.6
1.485956
1.485955
输出变量: 解向量元素
(4)具体算例及求解结果:
例:用列选主元法求解下列线性方程组(课本P65例3-3)
求解结果:
3、 分解法求解线性方程组
(1)算法原理:
求解线性方程组 时,当对 进行 分解,则等价于求解 ,这时可归结为利用递推计算相继求解两个三角形(系数矩阵为三角矩阵)方程组,用顺代,由
求出 ,再利用回带,由 求出 。
例:根据给定的函数 的实例数据表,试用最小二乘法求二次拟合多项式。(课本P186习题3)
求解结果:
6、变步长梯形求积分
(1)算法原理:
设将积分区间 分成 等份,即有 个子区间,分点 ,其中步长
对于子区间 ,利用体型求其积分近似值
对于子区间 有
对于子区间 再取其中点
作新节点,此时区间数增加了一倍为 ,
0.2
1.184097
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
数值计算方法上机实验报告班级:农电0801 姓名:杨 昆学号:200801090122 算法原理以及程序框图①牛顿法求非线性方程● 算法原理:对于非线性方程()0f x =,若已知跟*x 的一个近似值k x ,将()f x 在k x 处展开成一阶泰勒公式:"'2()()()()()()2!k k k k f f x f x f x x x x x ξ=+-+-忽略高次项,有'()()()()k k k f x f x f x x x ≈+-右端是直线方程,用这个直线方程来近似非线性方程()f x 。
将非线性方程()0f x =的根*x 代入*()0f x =,即'*()()()0k k k f x f x x x +-≈解出*'()()k k k f x x x f x ≈-将右端取为1k x +,则1k x +是比k x 更接近于*x 的近似值,即1'()()k k k k f x x x f x +=-这就是牛顿迭代公式。
● 程序框图:●具体算例及求解结果:导出计算(0)c c>的牛顿迭代公式,并计算115。
迭代结果:10.75000010.72383710.723805●输入变量、输出变量说明:输入变量:0x迭代初值,ε迭代精度,N迭代最大次数输出变量:k当前迭代次数,x当前迭代值1②列主元消去法求解线性方程组●算法原理:高斯消去法是利用现行方程组初等变换中的一种变换,即用一个不为零的数乘一个方程后加只另一个方程,使方程组变成同解的上三角方程组,然后再自下而上对上三角方程组求解。
列选住院是当高斯消元到第k 步时,从k 列的kk a 以下(包括kk a )的各元素中选出绝对值最大的,然后通过行交换将其交换到kk a 的位置上。
交换系数矩阵中的两行(包括常数项),只相当于两个方程的位置交换了,因此,列选主元不影响求解的结果。
● 程序框图:● 具体算例及求解结果:用列选主元法求解下列线性方程组1231231232315410030.12x x x x x x x x x ++=⎧⎪++=⎨⎪-+=⎩ 求解结果:1.2000002.000000 -1.400000● 输入变量、输出变量说明:输入变量:ij a 系数矩阵元素,i b 常向量元素 输出变量:12,,n b b b 解向量元素③LU 分解求解线性方程组● 算法原理:求解线性方程组Ax b =时,当对A 进行杜里特尔分解,则等价于求解LUx b =,这时可归结为利用递推计算相继求解两个三角形(系数矩阵为三角矩阵)方程组,用顺代,由Ly b =求出y ,再利用会带,由Ux y =求出x 。
● 程序框图:● 具体算例及求解结果:用杜里特尔分解法求解方程组123223347712457x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦ 求解结果:2.000000 -2.000000 1.000000● 输入变量、输出变量说明:输入变量:ij a 系数矩阵元素,i b 常向量元素 输出变量:12,,n b b b 解向量元素④拉格朗日插值多项式● 算法原理:首先构造基函数0()nik i k ii k x x l x x x =≠-=-∏,可以证明基函数满足下列条件:()1k i i k l x i k≠⎧=⎨=⎩,对于给定1n +个节点,n 次拉格朗日插值多项式由下式给出:0()()nk k k L x l x y ==∑其中0()nik i k ii kx x l x x x =≠-=-∏由于()k l x 是一个关于x 的n 次多项式,所以()L x 为关于x 的不高于n 次的代数多项式。
当i x x =时,()i i L x y =,满足插值条件。
● 程序框图:● 具体算例及求解结果:已知()sin f x x = 的值如下表所示。
()sin f x x =的值x0 6π 4π 3π 2π sin x1222321试用四次Lagrange 多项式计算sin 12π的估计值。
求解结果: 0.258588● 输入变量、输出变量说明:输入变量:(,)i i x y 插值节点输出变量:y 插值所得到被插函数在插值点的近似值⑤最小二乘法曲线拟合● 算法原理:对于给定的一组数据(,)i i x y ,i =1,2…,m ,寻求做n 次多项式0nk k k y a x ==∑使性能指标2011(,,,)()mnk n i k i i k J a a a y a x ===-∑∑ 为最小。
由于性能指标J 可以被看做关于k a ,k =0,1,…,n 的多元函数,故上述拟合多项式的构造问题可转化为多元函数的极值问题。
令0kJa ∂=∂ 从而有正则方程组203121122n i iii n i i i i iin n n n n n i i i i ii m x xx y a x y a x x xxa x y x x xx +++⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦∑∑∑∑∑∑∑∑∑∑∑∑∑∑求解即得多项式系数。
● 程序框图:● 具体算例及求解结果:给定函数()y f x =的实例数据表:i x 1 2 3 4 6 7 8 i y2367532试用最小二乘法求二次拟合多项式 求解结果: -1.318171 3.431811 -0.38636321.318171 3.4318110.386363y x x =-+-● 输入变量、输出变量说明:输入变量:(,)i i x y 已知数据点 输出变量:i a 拟合多项式的系数⑥改进欧拉法求解常微分方程初值问题● 算法原理:当h 取值较小时,让梯形法的迭代公式只迭代一次就结束。
这样先用欧拉公式求得一个初步近似值(0)1n y +,称之为预报值,预报值的精度不高,用它替代梯形法右端的1n y +,再直接计算得出1n y +,并称之为校正值,这时得到预报-校正公式。
将预报-校正公式(0)1(0)111(,)(,)(,)2n n n n n n n n n n y y hf x y h y y f x y f x y ++++⎧=+⎪⎨⎡⎤=++⎪⎣⎦⎩ 称为改进欧拉公式。
● 程序框图:● 具体算例及求解结果:求解初值问题'2,01(0)1x y y x y y ⎧=-<<⎪⎨⎪=⎩求解结果:n xn y()n y xn xn y()n y x0.1 1.095909 1.095909 0.6 1.485956 1.485955 0.2 1.184097 1.184097 0.7 1.562514 1.552514 0.3 1.266201 1.266201 0.8 1.616475 1.616474 0.4 1.343360 1.343360 0.9 1.678320 1.678166 0.5 1.416402 1.416402 1.0 1.737867 1.737867● 输入变量、输出变量说明:输入变量:00(,)x y 处置点,h 区间长度,N 计算次数 输出变量:11(,)x y 初值问题的数值解法结果⑦四阶龙格-库塔法求解常微分方程初值问题算法原理:用区间[]1,k k x x +内四个不同点上的函数值的线性组合就得到四阶龙格-库塔法。
四阶龙格-库塔法1112233441211113221122243311322333()(,)(,)(,)(,)n n n n n n n n n n y y h k k k k k f x y k f x h y k h k f x h y k h k h k f x h y k h k h k h ωωωωλμλμμλμμμ+=++++⎧⎪=⎪⎪=++⎨⎪=+++⎪⎪=++++⎩ 其中,,,ωλμ均为待定系数。
类似于前面的讨论,把234,,k k k 分别在n x 点展开成h 的幂级数,代入1n y +并进行花间,然后与1()n y x +在n x 点上的泰勒展开式比较,使其两式比较,使其两式右端直到4h 的系数相等,经过复杂的数学演算可得到关于,,ωλμ的一组特解121122213132333142312011613λλμμμμμλμωωωω⎧⎪====⎪⎪===⎪⎪==⎨⎪⎪==⎪⎪==⎪⎩从而得到下列常用的经典公式11234121123122413(22)6(,)(,)2(,)2(,)n n n n nn nn n n hy y k k k k k f x y h k f x y k h k f x y k k f x y hk ++++⎧⎪=++++⎪⎪=⎪⎪=+⎨⎪⎪=+⎪⎪⎪=+⎩ 经典的龙格-库塔法每一步需要4次计算函数值(,)f x y ,它具有四阶精度,即局部截断误差是5()O h 。
● 程序框图:● 具体算例及求解结果:求解初值问题'2,01(0)1x y y x y y ⎧=-<<⎪⎨⎪=⎩求解结果:n xn y()n y x0.2 1.183229 1.183229 0.4 1.341667 1.341667 0.6 1.483281 1.483281 0.8 1.612514 1.612514 1.01.732142 1.732142● 输入变量、输出变量说明:输入变量:00(,)x y 处置点,h 区间长度,N 计算次数数值计算方法上机实验报告11输出变量:11(,)x y 初值问题的数值解法结果上机体验与收获首先,是程序的编制过程中,我体会到了将一个抽象的数值计算方法变为一个具体计算机程序的不易。
这个过程,让我对各种数值计算方法的数学公式有了更进一步的认识,因为要编制出计算机程序,就必须清楚数值计算方法的具体计算步骤。
并且这次我还补充了自己的C 语言知识,学会了C 语言中文件的操作,这一部分内容在原先的C 语言课程中式没有讲述过的。
在各个程序的调试过程中,我又重新复习C 语言的很多知识,比如数组元素的编号以及数组行列的编号都是从0开始的等。
数值计算方法编程过程中涉及到了许多循环结构,这很好的锻炼了我的数学思维。
其次,我要说的事此次试验报告的写作,运用了Mathtype 在Word 里面编辑了如此多的数学公式,让我对运用Mathtype 进行数学公式编辑以及Word 排版都有了新的认识。
参考文献[1]谭浩强,C 程序设计,北京,清华大学出版社,2005[2]马东升,雷勇军,数值计算方法,机械工业出版社,2006。