MATLAB上机作业(数值分析方法(上理))
数值分析(王能超)matlab上机作业
引论习题题目:用二分法求方程x^3-x-1=0在[1,2]内的近似根,要求误差不超过10-3程序:function [root,n]=dichotomy(a,b,eps)% dichotomy:二分法求根函数% f:带求解方程;[a,b]求解区间;n为二分次数if nargin<2 disp('输入变量不够');return;else if nargin>3 disp('输入变量过多');return;else nargin==2eps=1.0e-3;endendif a>b disp('区间输入错误')returnendc=(a+b)/2;if f(c)==0y=cn=1returnendn=0;while abs(b-a)>epsif f(a)*f(c)>0a=c;else b=c;endc=(a+b)/2;n=n+1;endroot=c;end运算结果:定义函数f如下:function [yy] = f(x)%定义了测试函数% x为输入变量yy=x^3-x-1;end在commond window中运行结果如下:>> [root,n]=dichotomy(1,2,1.0e-3)root =1.32471n =10习题一题目一给出概率积分2xxy e dx-=的数据表i 0 1 2 3xi 0.46 0.47 0.48 0.49yi 0.4846555 0.4937452 0.5027498 0.5116683 用二次插值计算,试问:(1)当x=0.472时该积分值等于多少?(2)当x为何值时积分值等于0.5?程序:function [ y,t ] = interpolation( X,Y,x )%拉格朗日插值函数%X为自变量数组,Y为函数值数组,x为插值点,t为插值基函数n=length(X);y=0;t=zeros(1,n);for i=1:nt(i)=1;for j=1:nif i==j continueendt(i)=t(i)*(x-X(j))/(X(i)-X(j));endy=y+t(i)*Y(i);endend运算结果:在commond window中运行结果如下:(1)当x=0.472时该积分值等于多少?>> clear>> X=[0.46 0.47 0.48 0.49];>> Y=[0.4846555 0.4937452 0.5027498 0.5116683];>> x=0.472;>> [ y,t ] = interpolation( X,Y,x )y =0.4956t =-0.0480 0.8640 0.2160 -0.0320(2)当x为何值时积分值等于0.5?>> clear>> X=[0.46 0.47 0.48 0.49];>> Y=[0.4846555 0.4937452 0.5027498 0.5116683];>> y=0.5;>> [ x,t ] = interpolation( Y,X,y )x =0.4769t =-0.0452 0.3356 0.7707 -0.0611题目二构造适合下列数据表的三次插值样条函数S(x):x -1 0 1 3 y -1 1 3 5 y' 6 1程序:function [ y,m ] =spline( X,Y,x,m1,mn )%样条插值函数%X,Y为输入的自变量、函数值数组;%x为插值点(数组),y为插值结果(数组),m为各插值节点的一阶导数值%m1为X(1)的一阶导数;mn为X(n)的一阶导数n=length(X);alpha=zeros(n-2,1);beta=zeros(n-2,1);for i=1:(n-2)alpha(i)=(X(i+1)-X(i))/(X(i+1)-X(i)+X(i+2)-X(i+1));beta(i)=3*((1-alpha(i))*(Y(i+1)-Y(i))/(X(i+1)-X(i))+alpha(i)*(Y(i+2)-Y(i+1))/(X(i+2)-X(i+1) ));endm=zeros(n,1);m(1)=m1;m(n)=mn;A=zeros((n-2),(n-2));B=zeros((n-2),1);for j=1:(n-2)A(j,j)=2;endfor k=2:(n-2)A(k,(k-1))=1-alpha(k);endfor l=1:(n-3)A(l,(l+1))=alpha(l);endB(1)=beta(1)-(1-alpha(1))*m1;B(n-2)=beta(n-2)-alpha(n-2)*mn;for p=2:(n-3)B(p)=beta(p);endm(2:(n-1))=A\B;s=length(x);y=zeros(1,s);for t=1:sr=1;for q=1:nif x(t)>=X(q)&&x(t)<=X(q+1)break;endr=r+1;endxx=(x(t)-X(r))/(X(r+1)-X(r));y(t)=(xx-1)^2*(2*xx+1)*Y(r)+xx^2*(-2*xx+3)*Y(r+1)+(X(r+1)-X(r))*xx*(xx-1)^2*m(r)+( X(r+1)-X(r))*xx^2*(xx-1)*m(r+1);endend运算结果:在commond window中运行结果如下:>> clear>> X=[-1 0 1 3];>> Y=[-1 1 3 5];>> m1=6;mn=1;>> x=-1:0.1:3;>> [ y,m ] =spline( X,Y,x,m1,mn );plot(x,y)求得结果如下(整理之后):其中红色数据表示插值节点将以上结果绘制成图,如下所示:习题二题目编写一通用型复化辛普生公式,能够对任意长度的等间距离散数据进行积分运算。
基于MATLAB的数值分析编程上机作业1
%引用函数:
% holder2;示例[p,u]=holder2(x);
%使用举例:
% H=holderk(x,k)
%Define variables:
% x-输入的n维向量;
% n-n维向量x的维数;
% p-Householder初等变换阵的系数ρ;
H(k:n,k:n)=eye(n-k+1)-p\u*u';%计算H(k:n,k:n)=I-p\u*u';
3、计算实例:
>> x=[2,0,2,1]'
x =
2
0
2
1
>> H=holderk(x,3)
H =
1.0000 0 0 0
0 1.0000 0 0
0 0 -0.8944 -0.4472
0 0 -0.4472 0.8944
%HOLDERK给定向量x≠0,数k,计算初等反射阵Hk,使HkX=Y,其中Y的第k+1项到最后项全为零;
%程序功能:函数holderk给定向量x≠0,数k,计算初等反射阵Hk,使HkX=Y,%程序功能:函数holder2给定向量x≠0,计算Householder初等变换阵的p,u;
%输入:n维向量x,数k;
>> H*x
ans =
2.0000
0
-2.2361
0
二、用Householder变换法求矩阵A的正交分解A=QR。
1、基本原理:
任一实列满秩的m×n矩阵A,可以分解成两个矩阵的乘积,即A=QR,其中Q是具有法正交列向量的m×n矩阵,R是非奇异的n阶上三角阵。
数值分析作业MATLAB
1.用二分法解方程 x-lnx=2 在区间【2 ,4】内的根方法: 二分法算法:f=inline('x-2-log(x)');a=2;b=4;er=b-a; ya=f(a);er0=.00001;while er>er0x0=.5*(a+b);y0=f(x0);if ya*y0<0b=x0;elsea=x0;ya=y0;enddisp([a,b]);er=b-a;k=k+1;end求解结果:>> answer13 43.0000 3.50003.0000 3.25003.1250 3.25003.1250 3.18753.1250 3.15633.1406 3.15633.1406 3.14843.1445 3.1484 3.1445 3.1465 3.1455 3.1465 3.1460 3.1465 3.1460 3.1462 3.1461 3.1462 3.1462 3.14623.1462 3.1462 3.1462 3.1462 3.1462 3.1462 最终结果为: 3.14622.试编写MATLAB 函数实现Newton 插值,要求能输出插值多项式。
对函数141)(2+=x x f 在区间[-5,5]上实现10次多项式插值。
Matlab 程序代码如下:%此函数实现y=1/(1+4*x^2)的n 次Newton 插值,n 由调用函数时指定 %函数输出为插值结果的系数向量(行向量)和插值多项式 算法:function [t y]=func5(n) x0=linspace(-5,5,n+1)'; y0=1./(1.+4.*x0.^2); b=zeros(1,n+1); for i=1:n+1 s=0; for j=1:i t=1; for k=1:iif k~=jt=(x0(j)-x0(k))*t;end;end;s=s+y0(j)/t;end;b(i)=s;end;t=linspace(0,0,n+1);for i=1:ns=linspace(0,0,n+1);s(n+1-i:n+1)=b(i+1).*poly(x0(1:i));t=t+s;end;t(n+1)=t(n+1)+b(1);y=poly2sym(t);10次插值运行结果:[b Y]=func5(10)b =Columns 1 through 4-0.0000 0.0000 0.0027 -0.0000Columns 5 through 8-0.0514 -0.0000 0.3920 -0.0000Columns 9 through 11-1.1433 0.0000 1.0000Y =- (7319042784910035*x^10)/147573952589676412928 + x^9/18446744073709551616 + (256*x^8)/93425 -x^7/1152921504606846976 -(28947735013693*x^6)/562949953421312 -(3*x^5)/72057594037927936 + (36624*x^4)/93425 -(5*x^3)/36028797018963968 -(5148893614132311*x^2)/4503599627370496 +(7*x)/36028797018963968 + 1b为插值多项式系数向量,Y为插值多项式。
数值分析上机题Matlab(东南大学)3
0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.38 0.40 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 0.60 0.62 0.64 0.66 0.68 0.70 0.72
152 139 128 119 110 103 96 90 85 80 76 72 68 65 62 59 56 53 51 49 47 45 43 41 39 38
========================================================================================================================
======================================================================================================================================================================== 习题 3_36 ======================================================================================================================================================================== Omega n x1 x2 x3 x4 x5 x6 x7 x8 x9
-0.71279 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281
数值分析作业-matlab上机作业
数值分析———Matlab上机作业学院:班级:老师:姓名:学号:第二章解线性方程组的直接解法第14题【解】1、编写一个追赶法的函数输入a,b,c,d输出结果x,均为数组形式function x=Zhuiganfa(a,b,c,d)%首先说明:追赶法是适用于三对角矩阵的线性方程组求解的方法,并不适用于其他类型矩阵。
%定义三对角矩阵A的各组成单元。
方程为Ax=d%b为A的对角线元素(1~n),a为-1对角线元素(2~n),c为+1对角线元素(1~n-1)。
% A=[2 -1 0 0% -1 3 -2 0% 0 -2 4 -3% 0 0 -3 5]% a=[-1 -2 -3];c=[-1 -2 -3];b=[2 3 4 5];d=[6 1 -2 1];n=length(b);u(1)=b(1);y(1)=d(1);for i=2:nl(i)=a(i-1)/u(i-1);%先求l(i)u(i)=b(i)-c(i-1)*l(i);%再求u(i)%A=LU,Ax=LUx=d,y=Ux,%Ly=d,由于L是下三角矩阵,对角线均为1,所以可求y(i)y(i)=d(i)-l(i)*y(i-1);endx(n)=y(n)/u(n);for i=(n-1):-1:1%Ux=y,由于U是上三角矩阵,所以可求x(i)x(i)=(y(i)-c(i)*x(i+1))/u(i);end2、输入已知参数>>a=[2 2 2 2 2 2 2];>>b=[2 5 5 5 5 5 5 5];>>c=[2 2 2 2 2 2 2];>>d=[220/27 0 0 0 0 0 0 0];3、按定义格式调用函数>>x=zhuiganfa(a,b,c,d)4、输出结果x=[8.147775166909105 -4.073701092835030 2.036477565178471 -1.017492820111148 0.507254485099400 -0.250643392637350 0.119353996493976 -0.047741598597591]第15题【解】1、编写一个程序生成题目条件生成线性方程组A x=b 的系数矩阵A 和右端项量b ,分别定义矩阵A 、B 、a 、b 分别表示系数矩阵,其中1(10.1;,1,2,...,)j ij i i a x x i i j n -==+=或1(,1,2,...,)1ij a i j n i j ==+-分别构成A 、B 对应右端项量分别a 、b 。
数值计算上机实习题目(matlab编程)
数值计算上机实习题目(matlab编程)非线性方程求根一、实验目的本次实验通过上机实习,了解迭代法求解非线性方程数值解的过程和步骤。
二、实验要求1、用迭代法求方程230x x e -=的根。
要求:确定迭代函数?(x),使得x=?(x),并求一根。
提示:构造迭代函数2ln(3)x ?=。
2、对上面的方程用牛顿迭代计算。
3、用割线法求方程3()310f x x x =--=在02x =附近的根。
误差限为410-,取012, 1.9x x ==。
三、实验内容1、(1)首先编写迭代函数,记为iterate.mfunction y=iterate(x)x1=g(x); % x 为初始值。
n=1;while(abs(x1-x)>=1.0e-6)&(n<=1000) % 迭代终止的原则。
x=x1;x1=g(x);n=n+1;endx1 %近似根n %迭代步数(2)后编制函数文件?(x),记为g.mfunction y=g(x)y=log(3*x.^2);(3)设初始值为0、3、-3、1000,观察初始值对求解的影响。
将结果记录在文档中。
>>iterate(0)>>iterate(3) 等等2、(1)首先编制牛顿迭代函数如下,记为newton.mfunction y=newton(x0)x1=x0-fc(x0)/df(x0); % 牛顿迭代格式n=1;while(abs(x1-x0)>=1.0e-6)&(n<=1000000) % 迭代终止的原则。
x0=x1;x1=x0-fc(x0)/df(x0);n=n+1;endx1 %近似根n %迭代步数(2)对题目中的方程编制函数文件,记为fc.mfunction y=fc(x)y=3*x.^2-exp(x)编制函数的导数文件,记为df.mfunction y=df(x)y=6*x-exp(x)(3)在MATLAB 命令窗计算,当设初始值为0时,newton(0);给定不同的初始值,观察用牛顿法求解时所需要的迭代步数,并与上面第一题的迭代步数比较。
Matlab第一次上机作业
输入: >>tic, n=9;[u,k]=xsj(n), toc,surf(u)
计算Байду номын сангаас果如下:
u=
0 0 0 0 0 0 0 0 0 0 0 0 0.0256 0.0413 0.0508 0.0560 0.0577 0.0560 0.0508 0.0413 0.0256 0 0 0.0413 0.0686 0.0859 0.0955 0.0986 0.0955 0.0859 0.0686 0.0413 0 0.0508 0.0859 0.1088 0.1216 0.1258 0.1216 0.1088 0.0859 0.0508 0 0 0.0560 0.0955 0.1216 0.1364 0.1412 0.1364 0.1216 0.0955 0.0560 0 0 0.0577 0.0986 0.1258 0.1412 0.1462 0.1412 0.1258 0.0986 0.0577 0 0 0.0560 0.0955 0.1216 0.1364 0.1412 0.1364 0.1216 0.0955 0.0560 0 0 0.0508 0.0859 0.1088 0.1216 0.1258 0.1216 0.1088 0.0859 0.0508 0 0 0.0413 0.0686 0.0859 0.0955 0.0986 0.0955 0.0859 0.0686 0.0413 0 0 0.0256 0.0413 0.0508 0.0560 0.0577 0.0560 0.0508 0.0413 0.0256 0 0 0 0 0 0 0 0 0 0 0 0 0
程序一: function [u,k]=xsbj(n) % xsbj:用块 Jacobi 迭代法求解线性方程组 A*u=f % u:方程组的解; k 迭代次数; n:非边界点数 % a:方程组系数矩阵 Aii 的下对角线元素;b:方程系数矩阵 Aii 的主对角线元素 % c:方程组系数矩阵 Aii 的上对角线元素;d:追赶法所求方程的右端向量 % e:允许误差界;er:迭代误差 f=2*1/(n+1)^(2)*ones(n+2,n+2); a=-1*ones(1,n);b=4*ones(1,n);c=-1*ones(1,n);u=zeros(n+2,n+2);e=0.000000001; for k=1:2000 er=0; ub=u; for j=2:n+1 d(1:n)=f(2:n+1,j)+ub(2:n+1,j-1)+ub(2:n+1,j+1); x=zg(a,b,c,d); u(2:n+1,j)=x'; er=er+norm(ub(:,j)-u(:,j),1); end if er/n^2<e,break; end end 程序二: function x=zg(a,b,c,d) % zg:用追赶法求解线性方程组 n=length(b); % LU 分解 u(1)=b(1); for k=2:n if u(k-1)==0,D=0,return; end l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*c(k-1); end % 追赶法求解之追过程,求解 Ly=d y(1)=d(1); for k=2:n y(k)=d(k)-l(k)*y(k-1); end % 追赶法求解之追过程,求解 Ux=y if u(n)==0,D=0,return;end x(n)=y(n)/u(n); for k=n-1:-1:1 x(k)=(y(k)-c(k)*x(k+1))/u(k); end 输入:
MATLAB上机内容及作业
MATLAB上机内容及作业无约束优化求解函数fminsearch和fminunc求解无约束非线性优化问题的函数包括fminsearch函数和fminunc函数。
fminsearch和fminunc的函数相同,但fminunc函数可以在最优解处获得目标函数的梯度和Hessian矩阵值。
无约束优化数学模型为:minf(x)x∈rn求解无约束非线性优化问题的步骤为:第一步:先写目标函数的m文件;第二步是在命令窗口中调用相应的优化函数。
1.Fminsearch函数调用格式为[x,fval]=fminsearch(@myfun,x0)输出参数的含义:x:返回最优解的设计变量值;fval:在最优设计变量值时,目标函数的最小值;Exitflag:返回算法的终止标志。
在下列情况下,>0表示算法收敛到最优解;=0表示算法已达到最大迭代次数;<0表示算法不收敛。
output:返回优化算法信息的一个数据结构,有以下信息:输出Iteration表示迭代次数output Algorithm表示算法output Funccount表示函数求值的次数输入参数的含义:@Myfun:目标函数的m文件,前面加“@”或表示为“Myfun”,可以任意命名;X0:调用优化函数时,需要先给设计变量赋值;2、fminunc函数的调用格式[x,fval]=fminunc(@myfun,x0)grad:返回目标函数在最优解处的梯度信息;hessian:返回目标函数在最优解处的hessian矩阵信息。
其余含义同上。
3、实例已知某一优化问题的数学模型为:Minf(x)=3x12+2x1x2+x22x∈ r n由MATLAB程序编写的代码为:第一步:首先编写目标函数的.m文件,并保存为examplefsearch.m文件(先单击file菜单,后点击new命令中的m―file,即可打开m文件编辑窗口进行代码的编辑,在英文状态下输入程序代码),代码为:函数f=examplefsearch(x)f=3*x(1)^2+2*x(1)*x(2)+x(2)^2;1.0e-0.08*-0.79140.2260%分别为x1和x2的最优点的值(近似为0)Fval=1.5722e-016%,对应于最佳优势的最佳目标函数值(约为0)4。
数值分析MATLAB上机实验
数值分析实习报告姓名:gestepoA学号:201*******班级:***班序言随着计算机技术的迅速发展,数值分析在工程技术领域中的应用越来越广泛,并且成为数学与计算机之间的桥梁。
要解决工程问题,往往需要处理很多数学模型,不仅要研究各种数学问题的数值解法,同时也要分析所用的数值解法在理论上的合理性,如解法所产生的误差能否满足精度要求:解法是否稳定、是否收敛及熟练的速度等。
而且还能减少大量的人工计算。
由于工程实际中所遇到的数学模型求解过程迭代次数很多,计算量很大,所以需要借助如MATLAB,C++,VB,JAVA的辅助软件来解决,得到一个满足误差限的解。
本文所计算题目,均采用MATLAB进行编程,MATLAB被称为第四代计算机语言,利用其丰富的函数资源,使编程人员从繁琐的程序代码中解放出来MATLAB最突出的特点就是简洁,它用更直观的、符合人们思维习惯的代码。
它具有以下优点:1友好的工作平台和编程环境。
MATLAB界面精致,人机交互性强,操作简单。
2简单易用的程序语言。
MATLAB是一个高级的矩阵/阵列语言,包含控制语言、函数、数据结构,具有输入、输出和面向对象编程特点。
用户可以在命令窗口中将输入语句与执行命令同步,也可以先编好一个较大的复杂的应用程序(M文件)后再一起运行。
3强大的科学计算机数据处理能力。
包含大量计算算法的集合,拥有600多个工程中要用到的数学运算函数。
4出色的图像处理功能,可以方便地输出二维图像,便于我们绘制函数图像。
目录1 第一题.................................................................. 4 1.1 实验目的ﻩ41.2 实验原理和方法 (4)1.3 实验结果........................................................... 51.3.1 最佳平方逼近法ﻩ51.3.2 拉格朗日插值法ﻩ71.3.3 对比ﻩ82 第二题.................................................................. 92.1实验目的9ﻩ2.2实验原理和方法ﻩ10102.3实验结果ﻩ2.3.1 第一问10ﻩ2.3.2 第二问 (11)2.3.3 第三问1ﻩ13第三题 (12)3.1实验目的 (12)123.2 实验原理和方法ﻩ3.3实验结果1ﻩ24 MATLAB程序 (14)1第一题某过程涉及两变量x和y,拟分别用插值多项式和多项式拟合给出其对应规律的4.65880.3719.64484.272170345.2795 43 7.48478.2392 ⑴请用次数分别为3,4,5,6的多项式拟合并给出最好近似结果f(x)。
《MATLAB与数值分析》第一次上机实验报告
电子科技大学电子工程学院标准实验报告(实验)课程名称MATLAB与数值分析学生姓名:李培睿学号:2013020904026指导教师:程建一、实验名称《MATLAB与数值分析》第一次上机实验二、实验目的1. 熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算操作。
(用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序)2. 熟练掌握算术符号操作和基本运算操作,包括矩阵合并、向量合并、符号转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。
(用.m文件编写进行符号因式分解和函数求反的程序)3. 掌握Matlab函数的编写规范。
4、掌握Matlab常用的绘图处理操作,包括:基本平面图、图形注释命令、三维曲线和面的填充、三维等高线等。
(用.m文件编写在一个图形窗口上绘制正弦和余弦函数的图形,并给出充分的图形注释)5. 熟练操作MATLAB软件平台,能利用M文件完成MATLAB的程序设计。
三、实验内容1. 编程实现以下数列的图像,用户能输入不同的初始值以及系数。
并以x,y为坐标显示图像x(n+1) = a*x(n)-b*(y(n)-x(n)^2);y(n+1) = b*x(n)+a*(y(n)-x(n)^2)2. 编程实现奥运5环图,允许用户输入环的直径。
3. 实现对输入任意长度向量元素的冒泡排序的升序排列。
不允许使用sort 函数。
四、实验数据及结果分析题目一:①在Editor窗口编写函数代码如下:并将编写的函数文件用“draw.m”储存在指定地址;②在Command窗口输入如下命令:③得到图形结果如下:题目二:①在Editor窗口编写函数代码如下:并将编写的函数文件用“circle.m”储存在指定地址;②再次在Editor窗口编写代码:并将编写的函数文件用“Olympic.m”储存在指定地址;③在Command窗口输入如下指令(半径可任意输入):④按回车执行,将在图形窗口获得五环旗:题目三:①在Editor窗口编写函数代码如下:并用.将编写的函数文件用“qipaofa.m”储存在指定地址;②在Command窗口输入一组乱序数值,则可以得到升序排序结果如下:五、总结及心得体会1.要熟悉MATLAB编译软件的使用方法,明白有关语法,语句的基本用法,才可以在编写程序的时候游刃有余,不至于寸步难行。
数值分析matlab上机报告
不动点迭代:运行结果:代码截图:牛顿迭代:运行结果:代码截图:1)运行结果:代码截图:a二分法运行结果:代码截图:牛顿法:运行结果:代码截图:弦位法:代码截图:第二次上机作业运行结果:代码截图:A欧拉方法运行结果:代码截图:B改进的欧拉方法运行结果:代码截图:C四阶龙格库塔方法:代码截图:第三次上机作业第6章:(7版)P368: 5(e), 8, 9Using1) Gaussian elimination2) Gaussian elimination with partial pivoting.3) Gaussian elimination with scaled partial pivoting and three-digit chopping arithmetic to solvethe following linear systems, and compare the approximations to the actual solution.1.19x1 +2.11x2 − 100x3 + x4 = 1.12,14.2x1 − 0.122x2 + 12.2x3 − x4 = 3.44,100x2 − 99.9x3 + x4 = 2.15,15.3x1 + 0.110x2 − 13.1x3 − x4 = 4.16.Actual solution [0.176, 0.0126,−0.0206,−1.18]解:1)高斯消元法运行结果代码截图:2)部分选主元高斯消元法运行结果:代码截图:3)具有比例因子部分选主元高斯消元法运行结果:代码截图:2.Solve Ax = b using the Crout factorization for tridiagonal systems. Let A be the 10×10 tridiagonal matrix given byaii = 2, ai,i+1 = ai,i−1 = −1, for each i = 2, . . . , 9and a11 = a10,10 = 2, a12 = a10,9 = −1.Let b be the ten-dimensional column vector given byb1 = b10 = 1 and bi = 0, for each i = 2, 3, . . . , 9.运行结果:代码截图:3.1)雅克比迭代运行结果:程序截图:2)GS迭代运行结果代码截图:3)SOR方法运行结果:代码截图:。
数值分析上机作业
第二次上机作业一. 任务:用MATLAB 语言编写连续函数最佳平方逼近的算法程序(函数式M 文件)。
并用此程序进行数值试验,写出计算实习报告。
二. 程序功能要求:在后面的附一leastp.m 的基础上进行修改,使其更加完善。
要求算法程序可以适应不同的具体函数,具有一定的通用性。
所编程序具有以下功能:1. 用Lengendre 多项式做基,并适合于构造任意次数的最佳平方逼近多项式。
可利用递推关系 0112()1,()()(21)()(1)()/2,3,.....n n n P x P x xP x n xP x n P x n n --===---⎡⎤⎣⎦=2. 被逼近函数f(x)不用内联函数构造,而改用M 文件建立数学函数。
这样,此程序可通过修改建立数学函数的M 文件以适用不同的被逼近函数(要学会用函数句柄)。
3. 要考虑一般的情况]1,1[],[)(+-≠∈b a x f 。
因此,程序中要有变量代换的功能。
4. 计算组合系数时,计算函数的积分采用变步长复化梯形求积法(见附三)。
5. 程序中应包括帮助文本和必要的注释语句。
另外,程序中也要有必要的反馈信息。
6. 程序输入:(1)待求的被逼近函数值的数据点0x (可以是一个数值或向量)(2)区间端点:a,b 。
7. 程序输出:(1)拟合系数:012,,,...,n c c c c(2)待求的被逼近函数值00001102200()()()()()n n s x c P x c P x c P x c P x =++++三:数值试验要求:1. 试验函数:()cos ,[0,4]f x x x x =∈+;也可自选其它的试验函数。
2. 用所编程序直接进行计算,检测程序的正确性,并理解算法。
3. 分别求二次、三次、。
最佳平方逼近函数)x s (。
4. 分别作出逼近函数)x s (和被逼近函数)(x f 的曲线图进行比较。
(分别用绘图函数plot(0x ,s(0x ))和fplot(‘x cos x ’,[x 1 x 2,y 1,y 2]))四:计算实习报告要求:1.简述方法的基本原理,程序功能,使用说明。
数值分析作业MATLAB
数值分析作业MATLAB数值分析是研究用数值方法解决数学问题的一门学科。
它的主要目标是通过计算机编程解决数学问题,尤其是那些无法通过解析方法解决的问题。
MATLAB是一种常用的数值分析软件,它提供了丰富的数值计算函数和工具箱,能够方便地进行各种数值分析方法的实现和计算。
数值分析的研究内容很广泛,包括数值计算方法、数值逼近、数值微分和数值积分等。
在数值计算方法中,最常用的有数值解线性方程组、数值解非线性方程、数值积分、数值微分等。
例如,通过使用MATLAB的线性方程组求解函数或者工具箱中的线性代数函数,可以解决各种形式的线性方程组。
通过MATLAB的非线性方程求解函数,可以解决各种非线性方程的数值解。
而数值积分和数值微分则可以通过MATLAB的积分函数和微分函数来实现,实现对函数的积分和微分操作。
数值逼近是数值分析的重要内容之一,它研究的是如何用简单的函数逼近给定的复杂函数。
在MATLAB中,可以通过多项式逼近、三次样条、拉格朗日插值、最小二乘逼近等方法来实现数值逼近的计算。
例如,使用MATLAB的插值函数interp1可以实现一维函数的插值计算,使用MATLAB 的polyfit函数可以拟合一维数据集合的多项式曲线。
而对于二维函数和三维函数的逼近,可以使用MATLAB的interp2和interp3函数来实现。
数值微分和数值积分是数值分析中的基本操作之一、它们可以根据给定的函数计算函数的导数和积分。
在MATLAB中,使用diff函数可以计算一维函数的导数,使用trapz和quad函数可以计算一维函数的定积分和数值积分。
而对于二维函数和三维函数的微分和积分,可以使用MATLAB 的grad函数和integral2函数来实现。
此外,MATLAB还提供了很多其他的数学函数和工具,包括解微分方程、优化问题、曲线拟合和最小二乘等。
对于一些复杂的数学问题,可以通过使用MATLAB的符号计算工具箱来实现符号计算。
第四章数值分析上机作业
1.用两种方法解sin 5cos 20.31x x +=-,如何一次求出此方程的四个根。
解:方法一:在Matlab 命令窗口输入命令x=solve('sin(5*x)+cos(2*x)=-0.31')运行后输出结果:x =- 0.36685340479904406504913603551901 - 0.089925518994485866153169602644131*i 方法二:在Matlab 命令窗口输入命令E1=sym('sin(5*x)+cos(2*x)=-0.31');[x1]=solve(E1)运行后输出结果:x1 =- 0.36685340479904406504913603551901 - 0.089925518994485866153169602644131*i 2.用三种方法解方程11852912312x x x x -+-=。
解:方法一:在Matlab 命令窗口输入命令x=solve('9*x^11-12*x^8+x^5-3*x^2=12')运行后输出结果:x =1.1945355092112803561808169303487 - 0.25878939596021341127295614138703 + 0.98996423045277565907337492165509*i - 0.91081125361412082546165956220494 + 0.34557668006489020731472236114125*i - 0.6567748898900820562684890790666 - 0.94093805304063227438930031737768*i - 0.6567748898900820562684890790666 + 0.94093805304063227438930031737768*i 0.34695114229720670305614226506505 + 0.85428006847946699100354057810685*i - 0.25878939596021341127295614138703 - 0.98996423045277565907337492165509*i 0.88215664256156941185655405241917 + 0.47468310614563172153151897912015*i 0.34695114229720670305614226506505 - 0.85428006847946699100354057810685*i 0.88215664256156941185655405241917 - 0.47468310614563172153151897912015*i - 0.91081125361412082546165956220494 - 0.34557668006489020731472236114125*i 方法二:在Matlab 命令窗口输入命令fa=[9,0,0,-12,0,0,1,0,0,-3,0,-12];xk=roots(fa)运行后输出结果:xk =-0.9108 + 0.3456i-0.9108 - 0.3456i-0.6568 + 0.9409i-0.6568 - 0.9409i-0.2588 + 0.9900i-0.2588 - 0.9900i1.1945 + 0.0000i0.8822 + 0.4747i0.8822 - 0.4747i0.3470 + 0.8543i0.3470 - 0.8543i方法三:在Matlab命令窗口输入命令E1=sym('9*x^11-12*x^8+x^5-3*x^2-12=0'); [x1]=solve(E1)运行后输出结果:x1 =1.1945355092112803561808169303487 - 0.25878939596021341127295614138703 + 0.98996423045277565907337492165509*i - 0.91081125361412082546165956220494 + 0.34557668006489020731472236114125*i - 0.6567748898900820562684890790666 - 0.94093805304063227438930031737768*i - 0.6567748898900820562684890790666 + 0.94093805304063227438930031737768*i 0.34695114229720670305614226506505 + 0.85428006847946699100354057810685*i - 0.25878939596021341127295614138703 - 0.98996423045277565907337492165509*i 0.88215664256156941185655405241917 + 0.47468310614563172153151897912015*i 0.34695114229720670305614226506505 - 0.85428006847946699100354057810685*i 0.88215664256156941185655405241917 - 0.47468310614563172153151897912015*i - 0.91081125361412082546165956220494 - 0.34557668006489020731472236114125*3 用MATLAB 方法求函数11852()912312f x x x x x =-+--的导数'()f x 。
数值分析数学实验上机实验编程matlab源代码
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');ifstrcmp(ss,'y')x0=input('input initial value x0>>');k=0; elsebreakendendendk;%给出迭代次数x=x0;%给出解Gauss消元法源程序:cleara=input('输入系数阵:>>\n')b=input('输入列阵b:>>\n')n=length(b);A=[a b]x=zeros(n,1);%%%函数主体Yipusilong-1;for k=1:n-1;%%%是否进行主元选取if abs(A(k,k))<yipusilong;%事先给定的认为有必要选主元的小数yzhuyuan=1;elseyzhuyuan=0;endifyzhuyuan;%%%%选主元t=A(k,k);for r=k+1:n;if abs(A(r,k))>abs(t)p=r;else p=k;endend%%%交换元素if p~=k; for q=k:n+1;s=A(k,q);A(k,q)=A(p,q);A(p,q)=s;endendend%%%判断系数矩阵是否奇异或病态非常严重if abs(A(k,k))<yipusilongdisp(‘矩阵奇异,解可能不正确’)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;endendend%%%%求解xx(n)=A(n,n+1)/A(n,n);for k=n-1:-1:1;s=0;for r=k+1:n;s=s+A(k,r)*x(r);endt=(A(k,n+1)-s) x(k)=(A(k,n+1)-s)/A(k,k) end Lagrange插值源程序:n=input('将区间分为的等份数输入:\n');s=[-1+2/n*[0:n]];%%%给定的定点,Rf为给定的函数x=-1:0.01:1;f=0;for q=1:n+1;l=1;%求插值基函数for k=1:n+1;if k~=q;l=l.*(x-s(k))./(s(q)-s(k));elsel=l;endendf=f+Rf(s(q))*l;%求插值函数endplot(x,f,'r')%作出插值函数曲线grid onhold on分段线性插值源程序clearn=input('将区间分为的等份数输入:\n');s=[-1+2/n*[0:n]];%%%给定的定点,Rf为给定的函数m=0;hh=0.001;for x=-1:hh:1;ff=0;for k=1:n+1;%%%求插值基函数switch kcase 1if x<=s(2);l=(x-s(2))./(s(1)-s(2));elsel=0;endcase n+1if x>s(n);l=(x-s(n))./(s(n+1)-s(n)); elsel=0;endotherwiseif x>=s(k-1)&x<=s(k);l=(x-s(k-1))./(s(k)-s(k-1)); else if x>=s(k)&x<=s(k+1);l=(x-s(k+1))./(s(k)-s(k+1)); elsel=0;endendendff=ff+Rf(s(k))*l;%%求插值函数值endm=m+1;f(m)=ff;end%%%作出曲线x=-1:hh:1;plot(x,f,'r');grid onhold on三次样条插值源程序:(采用第一边界条件)clearn=input('将区间分为的等份数输入:\n');%%%插值区间a=-1;b=1;hh=0.001;%画图的步长s=[a+(b-a)/n*[0:n]];%%%给定的定点,Rf为给定的函数%%%%第一边界条件Rf"(-1),Rf"(1)v=5000*1/(1+25*a*a)^3-50/(1+25*a*a)^4;for k=1:n;%取出节点间距h(k)=s(k+1)-s(k);endfor k=1:n-1;%求出系数向量lamuda,miula(k)=h(k+1)/(h(k+1)+h(k));miu(k)=1-la(k);end%%%%赋值系数矩阵Afor k=1:n-1;for p=1:n-1;switch pcase kA(k,p)=2;case k-1A(k,p)=miu(p+1);case k+1A(k,p)=la(p-1);otherwiseA(k,p)=0;endendend%%%%求出d阵for k=1:n-1;switch kcase 1d(k)=6*f2c([s(k) s(k+1) s(k+2)])-miu(k)*v; case n-1d(k)=6*f2c([s(k) s(k+1) s(k+2)])-la(k)*v; otherwised(k)=6*f2c([s(k) s(k+1) s(k+2)]);endend %%%%求解M阵M=A\d';M=[v;M;v];%%%%m=0;f=0;for x=a:hh:b;if x==a;p=1;elsep=ceil((x-s(1))/((b-a)/n));endff1=0;ff2=0;ff3=0;ff4=0;m=m+1;ff1=1/h(p)*(s(p+1)-x)^3*M(p)/6;ff2=1/h(p)*(x-s(p))^3*M(p+1)/6;ff3=((Rf(s(p+1))-Rf(s(p)))/h(p)-h(p)*(M(p+1)-M(p))/6)*(x-s(p));ff4=Rf(s(p))-M(p)*h(p)*h(p)/6;f(m)=ff1+ff2+ff3+ff4 ;end%%%作出插值图形x=a:hh:b;plot(x,f,'k')hold on grid on 多项式最小二乘法源程序clear%%%给定测量数据点(s,f)s=[3 4 5 6 7 8 9];f=[2.01 2.98 3.50 5.02 5.47 6.02 7.05];%%%计算给定的数据点的数目n=length(f);%%%给定需要拟合的数据的最高次多项式的次数m=10;%%%程序主体for k=0:m;g=zeros(1,m+1);for j=0:m;t=0;for i=1:n;%计算内积(fai(si),fai(si))t=t+fai(s(i),j)*fai(s(i),k);endg(j+1)=t;endA(k+1,:)=g;%法方程的系数矩阵t=0;for i=1:n;%计算内积(f(si),fai(si))t=t+f(i)*fai(s(i),k);endb(k+1,1)=t;enda=A\b%求出多项式系数x=[s(1):0.01:s(n)]';y=0;fori=0:m;y=y+a(i+1)*fai(x,i);endplot(x,y)%作出拟合成的多项式的曲线grid onhold onplot(s,f,'rx') %在上图中标记给定的点表中,L10(x)为Lagrang e插值的10次多项式,S10(x),S40(x)分别代表n=10,40的三次样条插值函数,X10(x),X40(x)分别代表n=10,40的线性分段插值函数。
数值分析大作业(利用MATLAB软件)
实验报告课程名称:数值分析实验项目:曲线拟合/数值积分专业班级:姓名:学号:实验室号:实验组号:实验时间:20.10.24 批阅时间:指导教师:成绩:工业大学实验报告(适用计算机程序设计类)专业班级:学号:姓名:实验名称:曲线拟合与函数插值附件A 工业大学实验报告(适用计算机程序设计类)专业班级:学号:姓名:实验步骤或程序:附录一:1.利用二次,三次,四次多项式进行拟合:1.1 MATLAB代码如下:clear;clc;close allt=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24];y=[14 13 13 13 13 14 15 17 19 21 22 24 27 30 31 30 28 26 24 23 21 19 17 16 15];%输入数据hold on[p2 s2]=polyfit(t,y,2);%对于上面的数据进行2次多项式拟合,其中s2包括R(系数矩阵的QR分解的上三角阵),%df(自由度),normr(拟合误差平方和的算术平方根)。
y2=polyval(p2,t);%返回多项式拟合曲线在t处的值[p3 s3]=polyfit(t,y,3);y3=polyval(p3,t);[p4 s4]=polyfit(t,y,4);y4=polyval(p4,t);plot(t,y,'ro')%画图plot(t,y2,'g-')plot(t,y3,'m^-')plot(t,y4,'bs-')xlabel('t')ylabel('y')legend('原始数据','2次多项式拟合','3次多项式拟合','4次多项式拟合')1.2 二次,三次,四次多项式拟合的结果分别如下:(1)总的拟合结果在工作区的显示如下:(2)其次二次多项式拟合的结果为:(3)其中三次多项式拟合的结果:(4)其中四次多项式拟合的结果为:1.3 拟合的图像为:1.4 拟合的多项式为:根据工作区得出的数据列出最后的拟合多项式为:(1)y=7.416+2.594t-0.094t^2(2)y=12.251-0.102t+0.193t^2-0.008t^3(3)y=15.604-3.526t+0.866t^2-0.052t^3+0.0009t^42.形如2()()b t c y t ae--=的函数,其中,,a b c 为待定常数。
数值分析上机实验报告
一、实验目的通过本次上机实验,掌握数值分析中常用的算法,如二分法、牛顿法、不动点迭代法、弦截法等,并能够运用这些算法解决实际问题。
同时,提高编程能力,加深对数值分析理论知识的理解。
二、实验环境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、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
解得 x=0.8776 牛顿法求解:
建立第二个函数文件 function f2=f2(x) f2=15*x^2+12*x;
计算代码
x(1)=-1.7;jingdu=0.002; x(2)=x(1)-(f(x(1))/f2(x(1)));n=2;
while abs(x(n)-x(n-1))>jingdu x(n+1)=x(n)-(f(x(n))/f2(x(n))); n=n+1; end x
自定义计算文件
X = [0 pi/6 pi/4 pi/3 pi/2]; Y = [0 0.5 0.7071 0.8660 1]; x = linspace(0,pi,50); M = 1; [y, R] = lagrange(X, Y, x, M); y1 = sin(x); %画图 errorbar(x,y,R,'.k') hold on plot(X, Y, 'or', x, y, '*r', x, y1, '-b'); axis([-0.5,3.5,-1,1.2]);%把坐标边界加大一点,不然坐标点遮 住坐标轴了
end
A(i)=x^(i-1);
% 计算 G for j=1:n+1 for i=1:n+1 for k=1:m+1 B(k)=W(k)*subs(A(i),x,k)*subs(A(j),x,k); end G(i,j)=sum(B); end end G % 计算 d for j=1:n+1 for k=1:m+1 B(k)=W(k)*F(k)*subs(A(j),x,k); end d(j)=sum(B); end d=d' % 求出拟合曲线的系数 G=G^-1; C=G*d; S=0; for i=1:n+1 S=S+C(i)*x^(i-1); end % 画图 scatter(X,F,'*b'); hold on; ezplot(S,[0,6]); xlabel('x') ylabel('y') title('拟合曲线') end
end C = A(n, n); q1 = abs(q1*(z-X(n))); for k = (n-1):-1:1 C = conv(C, poly(X(k))); d = length(C); C(d) = C(d) + A(k,k);%在最后一维,也就是常数项加 上新的差商 end y(t) = polyval(C,z);%向量 y 是向量 x 处的插值,误差限 R,n 次牛顿插值多项式 L 及其系数向量 C end L = poly2sym(C); R(t) = M * q1 / c1;
上述代码插值部分和数值积分均参考自 CSDN 的几位知名博主, 大家可以自行查找。本文档旨在分享我的数值分析上机作业,不 做商用,故不与人产生任商业用途,转载请勿商用。想要系统 学习 matlab 软件,推荐飞天小苗苗的博客以及优酷她的教学视 频。
《数值计算方法》 上机报告
班级:机械*班 学号:171441**** 姓名:范涛
目 录 1. 非线性方程的数值解法 1.1 问题描述 1.2 基本原理和步骤 1.3 程序源代码 1.4 计算结果及分析 1.5 小结 2. 函数插值与曲线拟合 2.1 问题描述 2.2 基本原理和步骤 2.3 程序源代码 2.4 计算结果及分析 2.5 小结 3. 数值积分 3.1 问题描述 3.2 基本原理和步骤 3.3 程序源代码 3.4 计算结果及分析 3.5 小结
解得 x=0.8776 割线法求解: 计算代码
x(1)=-1.7;jingdu=0.002; x(2)=x(1)-(f(x(1))/f2(x(1)));n=2; while abs(x(n)-x(n-1))>jingdu x(n+1)=x(n)-(f(x(n))/(f(x(n))-f(x(n-1))))*(x(n)-x( n-1)); n=n+1; end x
解得 x=0.8776 1.4 计算结果分析 各种方法只要能够求解的情况下,多次求解就可以无限逼近所求解, 但是迭代法必须要求在一定范围内有极限, 而割线法要求一定范围内 必须为可切曲线,故有局限性 2.函数插值与曲线拟合 2.1 问题描述
用 Lagrange 插值来求 sinx 在某点的值,并估计其误差,已知 sin0° = 0, sin30° = 0.5, sin45° = 0.7071, sin60° = 0.8660, sin90° = 1.
解得 ans =18.849555921538759 复化三点式(辛普森公式)
function Sn = Sn(a,b,n) format long h = (b-a)/n; sum1 = 0; sum2 = 0; for i = 0:n-1 sum1 = sum1 + f(a+(i+1/2).*h); end for j = 1:n-1 sum2 = sum2 + f(a+j.*h); end Sn = h/6*(f(a)+4*sum1+2*sum2+f(b));
割线法:是利用牛顿迭代法的思想,在根的某个领域内,函数有直至 二阶的连续导数,并且不等于 0,则在领域内选取初值 x0,x1,迭代 均收敛。 1.3 程序源代码 建立被求解函数文件 function f=f(x) f=5*x^3+6*x^2-8; 二分法计算代码: a=-2;b=1;jingdu=0.002 while abs(a-b)>jingdu x=(a+b)/2;
c1 = c1 * j;
自定义计算文件
X = [0 pi/6 pi/4 pi/3 pi/2];%X 是 n+1 个节点(x_i,y_i)(i = 1,2, ... , n+1)横坐标向量,Y 是纵坐标向量 Y = [0 0.5 0.7071 0.8660 1]; x = linspace(0,pi,50); %x 是以向量形式输入的 m 个插值点,M
1.非线性方程的数值解法 1.1 令 f(x)=5*x^3+6*x^2-8=0,求根,x∈[-2,1],ε=0.002 1.2 基本原理和步骤 二分法:将[a,b]分为相等的 2 个小区间,计算小区间端点函数值, 找出新的有根区间。重复上述过程,直到找到满足精度要求的根。 迭代法:是取 x0 之后,在这个基础上,找到比 x0 更接近的方程的 跟,一步一步迭代,从而找到更接近方程根的近似跟。 牛顿法: 把非线性函数 f(x)=0 在 x0 处展开成泰勒级数取其线性 部分,作为非线性方程的近似方程, 则有 设 ,则其解为
legend('误差','样本点','拉格朗日插值估算','sin(x)'); title('拉格朗日插值法拟合函数');
牛顿插值法: 定义牛顿差值多项式函数
function[y,R,A,C,L] = newton(X,Y,x,M) n = length(X); m = length(x); for t = 1 : m z = x(t); A = zeros(n,n);%差商的矩阵 A A(:,1) = Y'; s = 0.0; p = 1.0; q1 = 1.0; c1 = 1.0; for j = 2 : n for i = j : n A(i,j) = (A(i,j-1) A(i-1,j-1))/(X(i)-X(i-j+1)); end q1 = abs(q1*(z-X(j-1)));
解得 ans =18.849555921538759 3.4 计算结果及分析 原函数需要多次积分能得到不定积分,应该是我的五次方对于我 的区间选择的影响较大,使得五次方之后最终结果前面的有效数 字接近,尝试把要求的函数降次积分之后返现结果误差较大,是 题目的问题,程序正确。 复化 Cotes 公式的代数精度最高, 复化 Simpson 其次, 复化梯形 公式的代数精度最低
2.2 基本原理和步骤 2.3 程序源代码 拉格朗日插值法: 拉格朗日插值函数文件
function [y, R] = lagrange(X, Y, x, M) n = length(X);m = length(x); for i = 1:m z = x(i); s = 0.0; for k = 1:n p = 1.0; q1 = 1.0; c1 = 1.0; for j = 1:n if j~=k p = p * (z - X(j)) / (X(k) - X(j)); end q1 = abs(q1 * (z - X(j))); c1 = c1 * j; end s = p * Y(k) + s; end y(i) = s; R(i) = M * q1 / c1; end
解得 ans =18.849555921538766 复化五点式(柯特斯公式)
function Cn = Cn(a,b,n) format long h = (b-a)/n;
sum1 = 0; sum2 = 0; for i = 0:n-1 sum1 = sum1 + 32*f(a+(i+1/4).*h)+12*f(a+(i+1/2).*h)+32*f(a+(i+3/4). *h); end for j = 1:n-1 sum2 = sum2 + 14*f(a+j.*h); end Cn = h/90*(7*f(a)+sum1+sum2+7*f(b));
命令行输入: X = [0 pi/6 pi/4 pi/3 pi/2];
F = [0 0.5 0.7071 0.8660 1]; W=[1 1 1 1 1]; S=mypolyfit(X,F,W,4,3)