数值分析MATLAB实验程序
MATLAB数值实验一(数据的插值运算及其应用完整版)
佛山科学技术学院实 验 报 告课程名称 数值分析 实验项目 插值法与数据拟合 专业班级 机械工程 姓 名 余红杰 学 号 10 指导教师 陈剑 成 绩 日 期 月 日一、实验目的1、学会Lagrange 插值、牛顿插值和三次样条插值等基本插值方法；2、讨论插值的Runge 现象3、学会Matlab 提供的插值函数的使用方法，会用这些函数解决实际问题。
二、实验原理1、拉格朗日插值多项式2、牛顿插值多项式3、三次样条插值 三、实验步骤1、用MATLAB 编写独立的拉格朗日插值多项式函数2、用MATLAB 编写独立的牛顿插值多项式函数3、用MATLAB 编写独立的三次样条函数（边界条件为第一、二种情形）4、已知函数在下列各点的值为：根据步骤1,2,3编好的程序，试分别用4次拉格朗日多项式4()L x 、牛顿插值多项式4()P x 以及三次样条函数()S x (自然边界条件)对数据进行插值，并用图给出 ｛(,),0.20.08,0,1,2,,10i i i x y x i i =+=｝，4()L x 、4()P x 和()S x 。
5、在区间[-1,1]上分别取10,20n =用两组等距节点对龙格函数21(),(11)125f x x x=-≤≤+作多项式插值，对不同n 值，分别画出插值函数及()f x 的图形。
6、下列数据点的插值可以得到平方根函数的近似，在区间[0,64]上作图。
（1）用这9个点作8次多项式插值8()L x 。
（2）用三次样条（第一边界条件）程序求()S x 。
7、对于给函数21()125f x x =+在区间[-1,1]上取10.2(0,1,,10)i x i i =-+=，试求3次曲线拟合，试画出拟合曲线并打印出方程，与第5题的结果比较。
四、实验过程与结果：1、Lagrange 插值多项式源代码：function ya=lag(x,y,xa) %x 所有已知插值点 %y 插值点对应函数值 %xa 所求点，自变量 %ya 所求点插值估计量 ya=0; mu=1; %初始化%循环方式求L 系数，并求和： for i = 1:length(y) for j = 1:length(x) if i ~= jmu = mu * (xa - x(j) ) / ( x(i) - x(j) ); else continue end endya = ya + y(i) * mu ; mu = 1; end2、Newton 源代码：function ya = newton(x,y,xa) %x 所有已知插值点 %y 插值点对应函数值 %xa 所求点，自变量 %ya 所求点插值估计量 %建立系数零矩阵D 及初始化：D = zeros(length(x)-1);ya = y(1);xi = 1;%求出矩阵D，该矩阵第一行为牛顿插值多项式系数：for i=1:(length(x)-1)D(i,1) = (y(i+1) -y(i))/(x(i+1) -x(i));endfor j=2:(length(x)-1)for i=1:(length(x)-j)D(i,j) = (D(i+1,j-1) - D(i,j-1)) / (x(i+j) - x(i)); endend%xi为单个多项式（x-x(1))(x-x(2))...的值for i=1:(length(x)-1)for j=1:ixi = xi*(xa - x(j));endya = ya + D(1,i)*xi;xi = 1;end3、三次样条插值多项式（1）（第一边界条件）源代码：function y=yt1(x0,y0,f_0,f_n,x) _____________(1)%第一类边界条件下三次样条插值；%xi 所求点；%yi 所求点函数值；%x 已知插值点；%y 已知插值点函数值；%f_0左端点一次导数值；%f_n右端点一次导数值；n = length(x0);z = length(y0);h = zeros(n-1,1);k=zeros(n-2,1);l=zeros(n-2,1);S=2*eye(n);for i=1:n-1h(i)= x0(i+1)-x0(i);endfor i=1:n-2k(i)= h(i+1)/(h(i+1)+h(i));l(i)= 1-k(i);end%对于第一种边界条件：k = [1;k]; _______________________(2)l = [l;1]; _______________________(3)%构建系数矩阵S：for i = 1:n-1S(i,i+1) = k(i);S(i+1,i) = l(i);end%建立均差表：F=zeros(n-1,2);for i = 1:n-1F(i,1) = (y0(i+1)-y0(i))/(x0(i+1)-x0(i));endD = zeros(n-2,1);for i = 1:n-2F(i,2) = (F(i+1,1)-F(i,1))/(x0(i+2)-x0(i));D(i,1) = 6 * F(i,2);end%构建函数D：d0 = 6*(F(1,2)-f_0)/h(1); ___________(4)dn = 6*(f_n-F(n-1,2))/h(n-1); ___________(5)D = [d0;D;dn]; ______________(6)m= S\D;%寻找x所在位置，并求出对应插值：for i = 1:length(x)for j = 1:n-1if (x(i)<=x0(j+1))&(x(i)>=x0(j))y(i) =( m(j)*(x0(j+1)-x(i))^3)/(6*h(j))+...(m(j+1)*(x(i)-x0(j))^3)/(6*h(j))+...(y0(j)-(m(j)*h(j)^2)/6)*(x0(j+1)-x(i))/h(j)+... (y0(j+1)-(m(j+1)*h(j)^2)/6)*(x(i)-x0(j))/h(j) ; break;else continue;endendend（2）（自然边界条件）源代码：仅仅需要对上面部分标注的位置做如下修改:__(1):function y=yt2(x0,y0,x)__(2):k=[0;k]__(3):l=[l;0]__(4)+(5):删除—(6):D=[0:D:0]4、——————————————PS:另建了一个f方程文件，后面有一题也有用到。
数值分析上机作业(MATLAB)
将系数矩阵 A 分解为：A=L+U+D
Ax=b
⇔ (D + L +U)x = b ⇔ Dx = −(L + U )x + b ⇔ x = −D −1(L + U )x + D −1b x(k +1) = −D −1 (L + U ) x(k ) + D −1b
输入 A，b 和初始向量 x
迭代矩阵 BJ , BG
否
ρ(B) < 1?
按雅各比方法进行迭代
否
|| x (k+1) − x(k) ||< ε ?
按高斯-塞德尔法进行迭代
否
|| x(k+1) − x (k ) ||< ε ?
输出迭代结果
图 1 雅各布和高斯-赛德尔算法程序流程图
1.2 问题求解
按图 1 所示的程序流程，用 MATLAB 编写程序代码，具体见附录 1。解上述三个问题 如下
16
-0.72723528355328
0.80813484897616
0.25249261987171
17
-0.72729617968010
0.80805513082418
0.25253982509100
18
-0.72726173942623
0.80809395746552
0.25251408253388
0.80756312717373
8
-0.72715363032573
0.80789064377799
9
-0.72718652854079
数值分析第四章外推法计算数值微分MATLAB计算实验报告
数值分析第四章外推法计算数值微分MATLAB计算实验报告数值分析MATLAB计算实验报告姓名班级学号⼀、实验名称⽤MATLAB编程实现数值微分的外推法计算。
⼆、实验⽬的1．掌握数值微分和定义和外推法的计算过程；2．了解数值微分外推法的计算⽅法并且编写出与其算法对应的MATLAB程序代码；3．体会利⽤MATLAB软件进⾏数值计算。
三、实验内容⽤外推法计算f(x)=x2e?x在x=0.5的导数。
四、算法描述1.命名函数。
2.如果输⼊未知数少于四个，默认精度10^-33.描述T表矩阵坐标4.依次赋值计算 T表第⼀列5.根据数值微分计算公式求出T表矩阵的值6.若达到精度则运算结束，若未达到循环计算7.输出T表，得出的值就是导数值五、实验结果六、实验结果分析此实验通过MATLAB实现外推法数值微分计算，得到相应的数据，⽅便对数据进⾏分析。
从结果可以看出，当步长h=0.025时⽤中点微分公式只有3位有效数字，外推⼀次达到5位有效数字，外推两次达到9位有效数字。
七、附录（程序）function g=waituifa(fname,x,h,e)if nargin<4,e=1e-3;end;i=1;j=1;G(1,1)=(feval(fname,x+h)-feval(fname,x-h))/(2*h);G(i+1,1)=(feval(fname,x+h/2)-feval(fname,x-h/2))/h;G(i+1,j+1)=(4^j*G(i+1,j)-G(i,j))/(4^j-1);while abs(G(i+1,i+1)-G(i+1,i))>ei=i+1;G(i+1,1)=(feval(fname,x+h/2^i)-feval(fname,x-h/2^i))/(2*h/2^i); for j=1:iG(i+1,j+1)=((4^j)*G(i+1,j)-G(i,j))/(4^j-1);endendGg=G(i+1,i+1);。
第十章MATLAB的数值分析
• 第一个问题可归结为“已知函数在x0,x1,
– …,xn处的值，求函数在区间[x0,xn]内其它点处的值”，这 种问题适宜用插值方法解决。 – 插值问题可描述为：已知函数在x0,x1,…,xn处的值 y0,y1,…,yn，求函数p(x)，使p(xi) = yi。
• 但对第二个问题不宜用插值方法,因为600米已超出所 给数据范围，用插值函数外推插值区间外的数据会 产生较大的误差。
– Q1=prctile(w,25); – Q3=prctile(w,75); – prctile( )函数实现计算样本的百分位数功能
分布形态的测定
• 只用集中趋势和离中趋势来表示所有数据，难免不 够准确。分析总体次数的分布形态有助于识别整个 总体的数量特征。总体的分布形态可以从两个角度 考虑，一是分布的对称程度，另一个是分布的高低。 前者的测定参数称为偏度或偏斜度，后者的测定参 数称为峰度。 • 峰度是掌握分布形态的另一指标，它能描述分布的 平缓或陡峭程度。如果峰度数值等于零，说明分布 为正态；若峰度数值大于零，说明分布呈陡峭状态； 若峰度数值小于零，说明分布形态趋于平缓。
– 解决第二个问题的常用方法是，根据地面到井下 500 处的 数据求出瓦斯浓度与地面到井下距离x之间的近似函数关 系f(x), 由f(x)求井下600米处的瓦斯浓度。
• 插值函数过已知点，拟合函数不一定过已知点。通 常, 插值主要用于求函数值，而拟合的主要目的是求 函数关系。当然，某些问题既可以用插值也可以用 拟合。
插值方法-概述
• 为什么需要插值？
(1) 函数关系y=f(x)没有明确的表达式
(2) y=f(x)表达式复杂，不便于研究和使用
-20 -15
沉陷量/mm 下沉方向为"+"
数值分析(hilbert矩阵)病态线性方程组的求解matlab程序
（Hilbert 矩阵）病态线性方程组的求解理论分析表明，数值求解病态线性方程组很困难。
考虑求解如下的线性方程组的求解Hx = b ，期中H 是Hilbert 矩阵，()ij n n Hh ，11ij h i j ，i ，j = 1,2，…，n 1.估计矩阵的2条件数和阶数的关系2.对不同的n ，取(1,1,,1)nx K ?，分别用Gauss 消去，Jacobi 迭代，Gauss-seidel 迭代，SOR 迭代和共轭梯度法求解，比较结果。
3.结合计算结果，试讨论病态线性方程组的求解。
第1小题：condition.m %第1小题程序t1=20;%阶数n=20x1=1:t1;y1=1:t1;for i=1:t1H=hilb(i);y1(i)=log(cond(H));endplot(x1,y1);xlabel('阶数n');ylabel('2-条件数的对数（log(cond(H)）');title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');t2=200;%阶数n=200x2=1:t2;y2=1:t2;for i=1:t2H=hilb(i);y2(i)=log(cond(H));endplot(x2,y2);xlabel('阶数n');ylabel('2-条件数的对数（log(cond(H)）');title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');画出Hilbert 矩阵2-条件数的对数和阶数的关系n=200时n=20时从图中可以看出，1）在n小于等于13之前，图像近似直线log（cond（H））~1.519n-1.8332）在n大于13之后，图像趋于平缓，并在一定范围内上下波动，同时随着n的增加稍有上升的趋势第2小题：solve.m%m第2小题主程序N=4000;xGauss=zeros(N,1);xJacobi=zeros(N,1);xnJ=zeros(N,1);xGS=zeros(N,1);xnGS=zeros(N,1);xSOR=zeros(N,1);xnSOR=zeros(N,1);xCG=zeros(N,1);xnCG=zeros(N,1);for n=1:N;x=ones(n,1);t=1.1;%初始值偏差x0=t*x;%迭代初始值e=1.0e-8;%给定的误差A=hilb(n);b=A*x;max=100000000000;%可能最大的迭代次数w=0.5;%SOR迭代的松弛因子G=Gauss(A,b);[J,nJ]=Jacobi(A,b,x0,e,max);[GS,nGS]=G_S(A,b,x0,e,max);[S_R,nS_R]=SOR(A,b,x0,e,max,w);[C_G,nC_G]=CG(A,b,x0,e,max);normG=norm(G'-x);xGauss(n)=normG;normJ=norm(J-x);nJ;xJacobi(n)=normJ;xnJ(n)=nJ;normGS=norm(GS-x);nGS;xGS(n)=normGS;xnGS(n)=nGS;normS_R=norm(S_R-x);nS_R;xSOR(n)=normS_R;xnSOR(n)=nS_R;normC_G=norm(C_G-x);nC_G;xCG(n)=normC_G;xnCG(n)=nC_G;endGauss.m%Gauss消去法function x=Gauss(A,b)n=length(b);l=zeros(n,n);x=zeros(1,n);%消去过程for i=1:n-1for j=i+1:nl(j,i)=A(j,i)/A(i,i);for k=i:nA(j,k)=A(j,k)-l(j,i)*A(i,k);endb(j)=b(j)-l(j,i)*b(i);endend%回代过程x(n)=b(n)/A(n,n);for i=n-1:-1:1c=A(i,:).*x;x(i)=(b(i)-sum(c(i+1:n)))/A(i,i);endJacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m 可能最大的迭代次数function [x,n]=Jacobi(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Jacobi迭代次数过多，迭代可能不收敛');break;endendG_S.m%Gauss-Seidel迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=G_S(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Gauss-Seidel迭代次数过多，迭代可能不收敛');break;endendSOR.m%SOR超松弛迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子function [x,n]=SOR(A,b,x0,e,m,w)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('SOR超松弛迭代次数过多，迭代可能不收敛');break;endendCG.m%CG共轭梯度法,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=CG(A,b,x0,e,m)r=b-A*x0;p=r;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=1;while norm(r1)>ebelta=(r1'*r1)/(r'*r);p=r1+belta*p;r=r1;x0=x;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=n+1;if n>mdisp('CG共轭梯度法迭代次数过多，迭代可能不收敛');break;endend。
数值分析中求解非线性方程的MATLAB求解程序
数值分析中求解非线性方程的MATLAB求解程序1. fzero函数：fzero函数是MATLAB中最常用的求解非线性方程的函数之一、它使用了割线法、二分法和反复均值法等多种迭代算法来求解方程。
使用fzero函数可以很方便地求解单变量非线性方程和非线性方程组。
例如，要求解方程f(x) = 0，可以使用以下语法：``````2. fsolve函数：fsolve函数是MATLAB中求解多维非线性方程组的函数。
它是基于牛顿法的迭代算法来求解方程组。
使用fsolve函数可以非常方便地求解非线性方程组。
例如，要求解方程组F(x) = 0，可以使用以下语法：``````3. root函数：root函数是MATLAB中求解非线性方程组的函数之一、它采用牛顿法或拟牛顿法来求解方程组。
使用root函数可以非常方便地求解非线性方程组。
例如，要求解方程组F(x) = 0，可以使用以下语法：``````4. vpasolve函数：vpasolve函数是MATLAB中求解符号方程的函数。
它使用符号计算的方法来求解方程，可以得到精确的解。
vpasolve函数可以求解多变量非线性方程组和含有符号参数的非线性方程。
例如，要求解方程组F(x) = 0，可以使用以下语法：```x = vpasolve(F(x) == 0, x)```vpasolve函数会返回方程组的一个精确解x。
5. fsolve和lsqnonlin结合：在MATLAB中，可以将求解非线性方程转化为求解最小二乘问题的形式。
可以使用fsolve函数或lsqnonlin函数来求解最小二乘问题。
例如，要求解方程f(x) = 0，可以将其转化为最小二乘问题g(x) = min，然后使用fsolve或lsqnonlin函数来求解。
具体使用方法可以参考MATLAB官方文档。
6. Newton-Raphson法手动实现：除了使用MATLAB中的函数来求解非线性方程，还可以手动实现Newton-Raphson法来求解。
数值分析MATLAB编程——数值积分法
数值分析MATLAB编程——数值积分法1、调用函数--f.Mfunction y=f(x)%------------------------------------------------------------函数1 y=sqrt(4-sin(x)*sin(x));%------------------------------------------------------------函数2 %y=sin(x)/x;%if x==0% y=0;%end%------------------------------------------------------------函数3 %y=exp(x)/(4+x*x);%------------------------------------------------------------函数4 %y=(log(1+x))/(1+x*x);2、复合梯形公式--tixing.M%复合梯形公式clear alla=input('请输入积分下限:');b=input('请输入积分上限:');n=input('区间n等分：');h=(b-a)/n;x=a:h:b;T=0;for k=1:n;T=0.5*h*(f(x(k))+f(x(k+1)))+T;endT=vpa(T,8)3、复合Simpson公式--simpson.M%复合Simpson公式clear alla=input('请输入积分下限:');b=input('请输入积分上限:');n=input('区间n等分：');h=(b-a)/n;x=a:h:b;S=0;for k=1:n;xx=(x(k)+x(k+1))/2;S=(1/6)*h*(f(x(k))+4*f(xx)+f(x(k+1)))+S;endS=vpa(S,8)4、Romberg算法--romberg.M%Romberg算法clear alla=input('请输入积分下限:');b=input('请输入积分上限:');n=input('区间n等分：');num=0:n;R=[num'];h=b-a;T=h*(f(a)+f(b))/2;t(1)=T;for i=2:n+1;u=h/2;H=0;x=a+u;while x<b;H=H+f(x);x=x+h;endt(i)=(T+h*H)/2;T=t(i);h=u;endR=[R,t'];for i=2:n+1for j=n+1:-1:1if j>=it(j)=(4^(i-1)*t(j)-t(j-1))/(4^(i-1)-1);elset(j)=0;endendR=[R,t'];endR=vpa(R,8)R(n,n)5、变步长算法（以复化梯形公式为例）--tixing2.M%复合梯形公式，确定最佳步长format longclear alla=input('请输入积分下限:');b=input('请输入积分上限:');eps=input('请输入误差：');k=1;T1=(b-a)*(f(a)+f(b))/2;T2=(T1+(b-a)*(f((a+b)/2)))/2; while abs((T1-T2)/3)>=epsM=0;n=2^k;h=(b-a)/n;T1=T2;x=a:h:b;for i=1:n;xx=(x(i)+x(i+1))/2;M=M+f(xx);endT2=(T1+h*M)/2;k=k+1;endT=vpa(T2,8)n=2^k。
数值分析matlab实验报告
数值分析matlab实验报告《数值分析MATLAB实验报告》摘要：本实验报告基于MATLAB软件进行了数值分析实验，通过对不同数学问题的数值计算和分析，验证了数值分析方法的有效性和准确性。
实验结果表明，MATLAB在数值分析领域具有较高的应用价值和实用性。
一、引言数值分析是一门研究利用计算机进行数值计算和分析的学科，其应用范围涵盖了数学、物理、工程等多个领域。
MATLAB是一种常用的数值计算软件，具有强大的数值分析功能，能够进行高效、准确的数值计算和分析，因此在科学研究和工程实践中得到了广泛的应用。
二、实验目的本实验旨在通过MATLAB软件对数值分析方法进行实验验证，探究其在不同数学问题上的应用效果和准确性，为数值分析方法的实际应用提供参考和指导。
三、实验内容1. 利用MATLAB进行方程求解实验在该实验中，利用MATLAB对给定的方程进行求解，比较数值解和解析解的差异，验证数值解的准确性和可靠性。
2. 利用MATLAB进行数值积分实验通过MATLAB对给定函数进行数值积分，比较数值积分结果和解析积分结果，验证数值积分的精度和稳定性。
3. 利用MATLAB进行常微分方程数值解实验通过MATLAB对给定的常微分方程进行数值解，比较数值解和解析解的差异，验证数值解的准确性和可靠性。
四、实验结果与分析通过对以上实验内容的实际操作和分析，得出以下结论：1. 在方程求解实验中，MATLAB给出的数值解与解析解基本吻合，验证了MATLAB在方程求解方面的高准确性和可靠性。
2. 在数值积分实验中，MATLAB给出的数值积分结果与解析积分结果基本吻合，验证了MATLAB在数值积分方面的高精度和稳定性。
3. 在常微分方程数值解实验中，MATLAB给出的数值解与解析解基本吻合，验证了MATLAB在常微分方程数值解方面的高准确性和可靠性。
五、结论与展望本实验通过MATLAB软件对数值分析方法进行了实验验证，得出了数值分析方法在不同数学问题上的高准确性和可靠性。
数值分析matlab完整版实验报告
《数值分析》报告运用Matlab求解非线性方程的根学院：专业：班级：姓名：学号：1. 目的掌握非线性方程求根的方法，并选取实例运用MATLAB 软件进行算法的实现，分别用牛顿法、弦截法和抛物线法求非线性方程的根。
2. 报告选题报告选取《数值分析（第四版）》290页习题7作为研究对象，即求3()310f x x x =--=在02x =附近的根。
根的准确值* 1.87938524...x =，要求结果准确到四位有效数字。
（1） 用牛顿法；（2） 用弦截法，取02x =，1 1.9x =； （3） 用抛物线法，取01x =，13x =，22x =。
3. 理论基础 (1) 牛顿迭代法牛顿迭代法是一种特殊的不动点迭代法,其计算公式为1(),0,1,2,...'()k k k k f x x x k f x +=-=其迭代函数为()()'()f x x x f x ϕ=-牛顿迭代法的收敛速度，当(*)0,'(*)0,''(*)0f x f x f x =≠≠时，容易证明,'(*)0f x ≠，''(*)''(*)0'(*)f x x f x ϕ=≠，牛顿迭代法是平方收敛的，且12''(*)lim2'(*)k k ke f x e f x +→∞=。
（2）弦截法将牛顿迭代法中的'()k f x 用()f x 在1k x -，k x 处的一阶差商来代替,即可得弦截法111()()()()k k k k k k k f x x x x x f x f x ++-=--- 。
（3）抛物线法弦截法可以理解为用过11(,()),(())k k k k x f x x f x ---两点的直线方程的根近似替()0f x =的根。
若已知()0f x =的三个近似根k x ，1k x -，2k x -用过1122(,()),(,()),(,())k k k k k k x f x x f x x f x ----的抛物线方程的根近似代替()0f x =的根，所得的迭代法称为抛物线法,也称密勒(Muller)法。
数值分析五个题目的C语言及MATLAB程序
1.7917595,1.9459101,2.079445,2.1972246,2.3025851}; //函数声明 void canshu(double *h,double *r,double *u,double *d); void wanju(double *u,double *r,double *d,double *m); double yths(double *m,double *h,double p); double ytdhs(double *m,double *h,double p); //主函数 void main() {
break; end end 实验二 C 程序 #include"stdio.h" #include"math.h" //已知量 double x[10]={1,2,3,4,5,6,7,8,9,10}; double fd1=1,fd2=0.1; double fx[10]={0,0.69314718,1.0986123,1.3862944,1.6094378,
} Matlab 程序 function[h,r,u,d]=canshu(x,fx,h,r,u,d,fd1,fd2) %UNTITLED10 Summary of this function goes here % Detailed explanation goes here for i=1:9
h(i)=x(i+1)-x(i); a(i)=(fx(i)-fx(i+1))/(x(i)-x(i+1)); j=i-1;
end end function [q] =ytdhs(m,h,p,fx,x) %UNTITLED14 Summary of this function goes here % Detailed explanation goes here for i=1:10
数值分析实验报告matlab
数值分析实验报告matlab数值分析实验报告引言：数值分析是一门研究利用计算机数值方法解决数学问题的学科，它在科学计算、工程设计、金融分析等领域具有重要的应用价值。
本实验报告旨在通过使用MATLAB软件，探索数值分析的基本原理和方法，并通过实际案例加深对数值分析的理解。
一、误差分析在数值计算中，误差是无法避免的。
误差分析是数值分析中的重要一环，它帮助我们了解数值计算的准确性和稳定性。
在实验中，我们通过计算机模拟了一个简单的数学问题，并分别计算了绝对误差和相对误差。
通过比较不同算法的误差大小，我们可以选择最适合的算法来解决实际问题。
二、插值与拟合插值和拟合是数值分析中常用的方法，它们可以通过已知的数据点来推导出未知数据点的近似值。
在本实验中，我们通过MATLAB的插值函数和拟合函数，分别进行了插值和拟合的实验。
通过比较不同插值和拟合方法的结果，我们可以选择最适合的方法来处理实际问题。
三、数值积分数值积分是数值分析中的重要内容，它可以用来计算曲线下的面积或函数的积分值。
在实验中，我们通过MATLAB的数值积分函数，对一些简单的函数进行了积分计算。
通过比较数值积分和解析积分的结果，我们可以评估数值积分的准确性和稳定性，并选择最适合的积分方法来解决实际问题。
四、常微分方程的数值解法常微分方程是数值分析中的重要内容，它可以用来描述许多自然现象和工程问题。
在实验中，我们通过MATLAB的常微分方程求解函数，对一些简单的微分方程进行了数值解法的计算。
通过比较数值解和解析解的结果，我们可以评估数值解法的准确性和稳定性，并选择最适合的数值解法来解决实际问题。
五、线性方程组的数值解法线性方程组是数值分析中的经典问题，它在科学计算和工程设计中广泛应用。
在实验中，我们通过MATLAB的线性方程组求解函数，对一些简单的线性方程组进行了数值解法的计算。
通过比较数值解和解析解的结果，我们可以评估数值解法的准确性和稳定性，并选择最适合的数值解法来解决实际问题。
数值分析matlab上机实验报告
数值分析matlab上机实验报告matlab软件实验报告数学上机课实验报告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) ); y(n+1) = b*x(n)+a*(y(n)-x(n) )2. 编程实现奥运5环图，允许用户输入环的直径。
3. 实现对输入任意长度向量元素的冒泡排序的升序排列。
不允许使用sort函数。
四、实验数据及结果分析题目一：①在Editor窗口编写函数代码如下：并将编写的函数文件用“draw.m”储存在指定地址；②在Command窗口输入如下命令：③得到图形结果如下：题目二：①在Editor窗口编写函数代码如下：并将编写的函数文件用“circle.m”储存在指定地址；②再次在Editor窗口编写代码：并将编写的函数文件用“Olympic.m”储存在指定地址；③在Command窗口输入如下指令（半径可任意输入）：④按回车执行，将在图形窗口获得五环旗：题目三：①在Editor窗口编写函数代码如下：并用.将编写的函数文件用“qipaofa.m”储存在指定地址；②在Command窗口输入一组乱序数值，则可以得到升序排序结果如下：五、总结及心得体会1. 要熟悉MATLAB编译软件的使用方法，明白有关语法，语句的基本用法，才可以在编写程序的时候游刃有余，不至于寸步难行。
matlab_examples(数值分析)
Experiments in Finding Root of Equation
计算方法课程组
华中科技大学数学与统计学院
1
方程求根 —二分法
一、实验目的
1) 熟悉Matlab编程 ; 2) 应用 Matlab实现二分法，牛顿迭代法等求根算法 ;
二、二分法
基本思想: 二分法通过不断搜索有根区间，最终收缩为一 点。算法简单、容易且保证算法收敛。
答案: x=1.3247
f(x)=-7.6580e-005
5
作
业
参看范例代码(.m文件) 独立完成如下编程内容:
6
2
方程求根 —二分法
function [xvect,xdif,fx,nit]=bisect(fun,a,b,toll,nmax)
% % % % % % % % % % 求根算法:二分法 [xvect,xdif,fx,nit]=bisect(fun,a,b,toll,nmax) fun 求根函数名 [a,b] 最初的有根区间的范围 toll 精度,默认为10e-5 nmax 最大迭代次数 xvect 返回所得根 xdif 返回缩小的根区间的长度 fx 返回函数值 nit 返回满足要求的迭代次数
bisect _main.m
x=1:0.01:2;
y=x.^3-x-1; plot(xห้องสมุดไป่ตู้y);hold on;
运行结果:
plot(x,zeros(size(x)),'r-.');
fun=inline('x^3-x-1');
[xvect,xdif,fx,nit]=bisect(fun,1,2,0.005,100);
disp((['
数值分析中求解非线性方程的MATLAB求解程序（6种）
数值分析中求解非线性方程的MATLAB求解程序（6种）数值分析中求解非线性方程的MATLAB求解程序（6种）1.求解不动点function [k,p,err,P]=fixpt(g,p0,tol,max1)%求解方程x=g(x) 的近似值，初始值为p0%迭代式为Pn+1=g(Pn)%迭代条件为：在迭代范围内满足|k|<1(根及附近且包含初值)k为斜率P(1)=p0;for k=2:max1P(k)=feval(g,P(k-1));err=abs(P(k)-P(k-1));relerr=err/(abs(P(k))+eps);p=P(k);if (err<tol)|(relerr<tol)< p="">break;endendif k==max1disp('超过了最长的迭代次数')endP=P';2.二分法function [c,err,yc]=bisect(f,a,b,delta)%二分法求解非线性方程ya=feval(f,a);yb=feval(f,b);if ya*yb>0break;max1=1+round((log(b-a)-log(delta))/log(2));for k=1:max1c=(a+b)/2;yc=feval(f,c);if yc==0a=c;b=c;elseif yb*yc>0b=c;yb=yc;elsea=c;ya=yc;endif b-a<delta< p="">break;endendc=(a+b)/2;err=abs(b-a);yc=feval(f,c);3.试值法function [c,err,yc]=regula(f,a,b,delta,epsilon,max1) %试值法求解非线性方程%f(a)和飞(b)异号ya=feval(f,a);yb=feval(f,b);if ya*yb>0disp('Note:f(a)*f(b)>0');for k=1:max1dx=yb*(b-a)/(yb-ya);c=b-dx;ac=c-a;yc=feval(f,c);if yc==0break;elseif yb*yc>0b=c;yb=yc;elsea=c;ya=yc;enddx=min(abs(dx),ac);if abs(dx)<delta|abs(yc)<epsilon< p="">break;endendc;err=abs(b-a)/2;yc=feval(f,c);4.求解非线性方程根的近似位置function R=approot(X,epsilon)%求解根近似位置%为了粗估算方程f(x)=0在区间[a,b]的根的位置，%使用等间隔采样点(xk,f(xk))和如下的评定准则：%f(xk-1)与f(xk)符号相反，%或者|f(xk)|足够小且曲线y=f(x)的斜率在%(xk,f(xk))附近改变符号。
数值分析matlab实验报告
数值分析matlab实验报告数值分析 Matlab 实验报告一、实验目的数值分析是研究各种数学问题数值解法的学科，Matlab 则是一款功能强大的科学计算软件。
本次实验旨在通过使用 Matlab 解决一系列数值分析问题，加深对数值分析方法的理解和应用能力，掌握数值计算中的误差分析、数值逼近、数值积分与数值微分等基本概念和方法，并培养运用计算机解决实际数学问题的能力。
二、实验内容（一）误差分析在数值计算中，误差是不可避免的。
通过对给定函数进行计算，分析截断误差和舍入误差的影响。
例如，计算函数＄f(x) ＝＼sin(x)＄在＄x ＝ 05$ 附近的值，比较不同精度下的结果差异。
（二）数值逼近1、多项式插值使用拉格朗日插值法和牛顿插值法对给定的数据点进行插值，得到拟合多项式，并分析其误差。
2、曲线拟合采用最小二乘法对给定的数据进行线性和非线性曲线拟合，如多项式曲线拟合和指数曲线拟合。
（三）数值积分1、牛顿柯特斯公式实现梯形公式、辛普森公式和柯特斯公式，计算给定函数在特定区间上的积分值，并分析误差。
2、高斯求积公式使用高斯勒让德求积公式计算积分，比较其精度与牛顿柯特斯公式的差异。
（四）数值微分利用差商公式计算函数的数值导数，分析步长对结果的影响，探讨如何选择合适的步长以提高精度。
三、实验步骤（一）误差分析1、定义函数｀compute_sin_error` 来计算不同精度下的正弦函数值和误差。
｀｀｀matlabfunction value, error ＝ compute_sin_error(x, precision)true_value ＝ sin(x)；computed_value ＝ vpa(sin(x)， precision)；error ＝ abs(true_value computed_value)；end｀｀｀2、在主程序中调用该函数，分别设置不同的精度进行计算和分析。
（二）数值逼近1、拉格朗日插值法｀｀｀matlabfunction L ＝ lagrange_interpolation(x, y, xi)n ＝ length(x)；L ＝ 0;for i ＝ 1:nli ＝ 1;for j ＝ 1:nif j ～＝ ili ＝ li （xi x(j)）／（x(i) x(j)）；endendL ＝ L ＋ y(i) li;endend｀｀｀2、牛顿插值法｀｀｀matlabfunction N ＝ newton_interpolation(x, y, xi)n ＝ length(x)；％计算差商表D ＝ zeros(n, n)；D(：， 1) ＝ y'；for j ＝ 2:nfor i ＝ j:nD(i, j) ＝（D(i, j 1) D(i 1, j 1)）／（x(i) x(i j ＋ 1)）；endend％计算插值结果N ＝ D(1, 1)；term ＝ 1;for i ＝ 2:nterm ＝ term （xi x(i 1)）；N ＝ N ＋ D(i, i) term;endend｀｀｀3、曲线拟合｀｀｀matlab％线性最小二乘拟合p ＝ polyfit(x, y, 1)；y_fit_linear ＝ polyval(p, x)；％多项式曲线拟合p ＝ polyfit(x, y, n)；％ n 为多项式的次数y_fit_poly ＝ polyval(p, x)；％指数曲线拟合p ＝ fit(x, y, ＇exp1'）；y_fit_exp ＝ p(x)；｀｀｀（三）数值积分1、梯形公式｀｀｀matlabfunction T ＝ trapezoidal_rule(f, a, b, n)h ＝（b a) ／ n;x ＝ a:h:b;y ＝ f(x)；T ＝ h （（y(1) ＋ y(end)）／ 2 ＋ sum(y(2:end 1)））；end｀｀｀2、辛普森公式｀｀｀matlabfunction S ＝ simpson_rule(f, a, b, n)if mod(n, 2) ～＝ 0error(＇n 必须为偶数'）；endh ＝（b a) ／ n;x ＝ a:h:b;y ＝ f(x)；S ＝ h ／ 3 （y(1) ＋ 4 sum(y(2:2:end 1)）＋ 2 sum(y(3:2:end 2)）＋ y(end)）；end｀｀｀3、柯特斯公式｀｀｀matlabfunction C ＝ cotes_rule(f, a, b, n)h ＝（b a) ／ n;x ＝ a:h:b;y ＝ f(x)；w ＝ 7, 32, 12, 32, 7 ／ 90;C ＝ h sum(w y)；end｀｀｀4、高斯勒让德求积公式｀｀｀matlabfunction G ＝ gauss_legendre_integration(f, a, b)x, w ＝ gauss_legendre(5)；％选择适当的节点数t ＝（b a) ／ 2 x ＋（a ＋ b) ／ 2;G ＝（b a) ／ 2 sum(w f(t)）；end｀｀｀（四）数值微分｀｀｀matlabfunction dydx ＝ numerical_derivative(f, x, h)dydx ＝（f(x ＋ h) f(x h)）／（2 h)；end｀｀｀四、实验结果与分析（一）误差分析通过不同精度的计算，发现随着精度的提高，误差逐渐减小，但计算时间也相应增加。
数值分析在生活中的应用举例及Matlab实现
一、最小二乘法，用MATLAB实现1. 数值实例下面给定的是乌鲁木齐最近1个月早晨7：00左右（新疆时间）的天气预报所得到的温度，按照数据找出任意次曲线拟合方程和它的图像。
下面用MATLAB编程对上述数据进行最小二乘拟合。
下面用MATLAB编程对上述数据进行最小二乘拟合2、程序代码x=[1:1:30];y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1];a1=polyfit(x,y,3) %三次多项式拟合%a2= polyfit(x,y,9) %九次多项式拟合%a3= polyfit(x,y,15) %十五次多项式拟合%b1=polyval(a1,x)b2=polyval(a2,x)b3=polyval(a3,x)r1= sum((y-b1).^2) %三次多项式误差平方和%r2= sum((y-b2).^2) %九次次多项式误差平方和%r3= sum((y-b3).^2) %十五次多项式误差平方和%plot(x,y,'*') %用*画出x,y图像%hold onplot(x,b1, 'r') %用红色线画出x,b1图像%hold onplot(x,b2, 'g') %用绿色线画出x,b2图像%hold onplot(x,b3, 'b:o') %用蓝色o线画出x,b3图像%3、数值结果不同次数多项式拟合误差平方和为：r1=67.6659r2=20.1060r3=3.7952r1、r2、r3分别表示三次、九次、十五次多项式误差平方和。
4、拟合曲线如下图二、 线性方程组的求解（ 高斯-塞德尔迭代算法 ）1、实例： 求解线性方程组（见书P233页）⎪⎪⎩⎪⎪⎨⎧=++=-+=+-3612363311420238321321321x x x x x x x x x 记A x=b, 其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=363320,,12361114238321b x A x x x任取初始值()()Tx0000=，进行迭代。
数值分析matlab数值试验
实验一：误差传播及算法稳定性实验1.21、试验程序：function charpt1_2% 误差传播及算法稳定性实验clc;clear all;promps={'请选择递推关系式，若选E1=1/e,En=1-nEn-1,请输入1，若选EN=0，En-1=(1-En)/n,请输入2：'};I=1;while Iresult=inputdlg(promps,'charpt1_2',1,{'1'});Nb=str2num(char(result));if ((Nb~=1)|(Nb~=2))I=0;endend%%%%%%%%%%%%%%%%%I=1;while Iresult=inputdlg('请输入递推步数 n>=1:','charpt1_2',1,{'10'});steps=str2num(char(result));if (steps>0)&(steps==fix(steps)) %% 如果steps大于0且为整数I=0;endend%%%%%%%%%%%%%%%%%result=inputdlg('请输入计算中所采用的有效数字位数n:','charpt1_2',1,{'5'});Sd=str2num(char(result));format long %% 设置显示精度result=zeros(1,steps); %% 存储计算结果err=result; %% 存储计算的绝对误差值func=result; %% 存储用quadl计算的近似值%%%%%%%%%%%%%%%%%%% 用quadl计算积分近似值for n=1:stepsfun=@(x) x.^n.*exp(x-1);func(n)= quadl(fun,0,1);end%%%%%%%%%%%%%%%%%%% 用自定义算法计算if(Nb==1)digits(Sd);result(1)=subs(vpa(1/exp(1)));for n=2:stepsresult(n)=subs(vpa(1-n*result(n-1)));enderr=abs(result-func);elseif(Nb==2)digits(Sd);result(steps)=0;for n=(steps-1):-1:1result(n)=subs(vpa((1-result(n+1))/(n+1)));enderr=abs(result-func);end%%%%%%%%%%%%%%%%%%% 输出结果数值及图像clf;disp('库函数计算值:');disp(sprintf('%e ',func));disp('递推值:');disp(sprintf('%e ',result));disp('误差值:');disp(sprintf('%e ',err));if(Nb==1)plot([1:steps],result,'-rs',[1:steps],func,':k*',[1:steps],err,'-.bo' );elseif(Nb==2)plot([steps:-1:1],result,'-rs',[steps:-1:1],func,':k*',[steps:-1:1],e rr,'-.bo');endxlabel('第n步');ylabel('计算值');legend('自定义算法结果','库函数计算结果','误差值');grid on2、试验结果：选择递推关系式1，递推步数为10，有效数字为5位，计算结果如下：库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-0011.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002递推值:3.678800e-001 2.642400e-001 2.072800e-001 1.708800e-001 1.456000e-0011.264000e-001 1.152000e-001 7.840000e-0022.944000e-001 -1.944000e+000误差值:5.588280e-007 1.117662e-006 3.352927e-006 1.341222e-0056.705713e-005 4.023702e-004 2.816427e-003 2.253226e-002 2.027877e-001 2.027877e+00012345678910第n 步计算值选择递推关系式2，递推步数为10，有效数字为5位，计算结果如下： 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678800e-001 2.642400e-001 2.072800e-001 1.708900e-001 1.455300e-001 1.267900e-001 1.125000e-001 1.000000e-001 1.000000e-001 0.000000e+000 误差值:5.588280e-007 1.117662e-006 3.352927e-006 3.412224e-006 2.942873e-006 1.237016e-005 1.164270e-004 9.322618e-004 8.387707e-003 8.387707e-002第n 步计算值选择递推关系式1，递推步数为10，有效数字为6位，计算结果如下： 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678790e-001 2.642420e-001 2.072740e-001 1.709040e-001 1.454800e-001 1.271200e-001 1.101600e-001 1.187200e-001 -6.848000e-002 1.684800e+000 误差值:4.411720e-007 8.823378e-007 2.647073e-006 1.058778e-0055.294287e-005 3.176298e-004 2.223573e-003 1.778774e-002 1.600923e-001 1.600923e+00012345678910第n 步计算值选择递推关系式2，递推步数为10，有效数字为6位，计算结果如下： 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678800e-001 2.642410e-001 2.072770e-001 1.708930e-001 1.455360e-001 1.267860e-001 1.125000e-001 1.000000e-001 1.000000e-001 0.000000e+000 误差值:5.588280e-007 1.176622e-007 3.529274e-007 4.122239e-007 3.057127e-006 1.637016e-005 1.164270e-004 9.322618e-004 8.387707e-003 8.387707e-002第n 步计算值3、结果分析：很明显第二种递推式结果要比第一种好，式1在第七步后有明显误差，而式2在第三步后基本与近似解一致。
数值分析课程设计报告(MATLAB版)
（2）取右端向量 b 的三位有效数字得 b [1.83 1.08 0.783]T ，求方程组的准确 解 X ，并与 X 的数据 [1 1 1]T 作比较 。说明矩阵的病态性。
算法及相应结果： （1）在 MATLAB 命令窗口里输入如下命令： >> H=[1 1/2 1/3;1/2 1/3 1/4;1/3 1/4 1/5]; b=[11/6 13/12 47/60]'; >> x=H\b 回车得到结果为： x = 1.0000 1.0000 1.0000 (2)紧接着在上题基础上继续输入如下命令： >> c=[1.83 1.08 0.783]'; x1=H\c 回车得到如下结果： x1 = 1.0800 0.5400 1.4400
问题分析：考虑由直线段（2 个点）产生第一个图形（5 个点）的过程，设 P 1 和 P5 分别为原始直线段的两个端点。现在需要在直线段的中间依次插入三个点 。显然， P2 位于 P P2 , P3 , P4 产生第一次迭代的图形（图 1-4） 1 点右端直线段的三分 之一处， P4 点绕 P2 旋转 60 度（逆时针方向）而得到的，故可以处理为向量 P2 P4 经正交变换而得到向量 P2 P3 ，形成算法如下： （1） P2 P 1 (P 5 P 1) / 3 ； （2） P4 P 1 2( P 5 P 1) / 3 ； （3） P3 P2 ( P4 P2 ) AT ； 在算法的第三步中，A 为正交矩阵。
运行结果： 0.0884 0.0580 0.0431 0.0343 0.0285 0.0243 0.0212 0.0188 0.0169 0.0154 0.0141 0.0130 0.0120 0.0112 0.0105 0.0099 0.0094 0.0087 0.0092 0.0042 （2）从 I 30 较粗略的估计值出发，我们不妨取 0.01. 源程序：
matlab数值分析实验报告
matlab数值分析实验报告Matlab数值分析实验报告引言数值分析是一门研究利用计算机进行数值计算和模拟的学科，它在科学计算、工程技术和金融等领域有着广泛的应用。
本次实验报告将介绍在Matlab环境下进行的数值分析实验，包括数值微分、数值积分和线性方程组求解等内容。
一、数值微分数值微分是通过数值方法计算函数的导数，常用的数值微分方法有前向差分、后向差分和中心差分。
在Matlab中，可以使用diff函数来计算函数的导数。
例如，对于函数f(x)=x^2，在Matlab中可以使用如下代码进行数值微分的计算：```matlabsyms x;f = x^2;df = diff(f, x);```二、数值积分数值积分是通过数值方法计算函数的定积分，常用的数值积分方法有梯形法则、辛普森法则和龙贝格积分法。
在Matlab中，可以使用trapz、quad和integral等函数来进行数值积分的计算。
例如，对于函数f(x)=sin(x)，可以使用如下代码进行数值积分的计算：```matlabx = linspace(0, pi, 100);y = sin(x);integral_value = trapz(x, y);```三、线性方程组求解线性方程组求解是数值分析中的重要问题，常用的求解方法有高斯消元法和LU 分解法。
在Matlab中，可以使用\操作符来求解线性方程组。
例如，对于线性方程组Ax=b，可以使用如下代码进行求解：```matlabA = [1, 2; 3, 4];b = [5; 6];x = A\b;```四、实验结果与分析在本次实验中，我们分别使用Matlab进行了数值微分、数值积分和线性方程组求解的计算。
通过实验结果可以发现，Matlab提供了丰富的数值计算函数和工具，能够方便地进行数值分析的计算和求解。
数值微分的计算结果与解析解相比较，可以发现数值微分的误差随着步长的减小而减小，但是当步长过小时，数值微分的误差会受到舍入误差的影响。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析上机实验报告院系:XXX学号:XXX姓名:XXXX目录一．Gass-Jordan消去法 (1)二．雅克比迭代法 (2)三．拉格朗日多项式插值 (4)四．埃尔米特插值法 (5)五．复化梯形公式 (7)六．龙贝格求积公式 (8)七．欧拉法 (9)八．隐式欧拉法 (10)一．Gass-Jordan消去法Gass-Jordan消去法求解线性方程组的基本思想是将系数矩阵化为对角矩阵，这样可以直接得到方程组的解x1,x2…x n，因此也称为无回代消去法。
例：用列主元Gass-Jordan消去法求解方程组：2x1+1x2+x3=55x1-1x2+1x3=81x1-3x2-4x3=-4编写Gau_Jor 函数来实现求解：function [x,flag]=Gau_Jor(A,b)% 求线性方程组的列主元Gauss－Jordan消去法% A为方程组的系数矩阵；% b为方程组的右端项；% x为方程组的解；% flag为指标向量，flag＝'failure'表示计算失败，flag='OK'表示计算成功。
[n,m]=size(A);nb=length(b);% 当方程组行与列的维数不相等时，停止计算，并输出出错信息if n~=merror('The rows and columns of matrix A must be equal!');return;end% 当方程组与右端项的维数不匹配时，停止计算，并输出出错信息if m~=nberror('The columns of A must be equal the length of b!');return;end% 开始计算，先赋初值flag='OK';x=zeros(n,1);for k =1:n% 选主元max1=0;for i=k:nif abs(A(i,k))>max1max1=abs(A(i,k));r=i;endendif max1<1e-10falg='failure';return;end% 交换两行if 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;end%消元计算b(k)=b(k)/A(k,k);for j=k+1:nA(k,j)=A(k,j)/A(k,k);endfor i=1:nif i~=kfor j =k+1:nA(i,j)=A(i,j)-A(i,k)*A(k,j);endb(i)=b(i)-A(i,k)*b(k);endendend% 输出xfor i=1:nx(i)=b(i);end在命令窗口输入：clearA=[2 1 2;5 -1 1;1 -3 -4];b=[5 8 -4];x=Gau_Jor(A,b)运行程序，输出如下：x =1.0000-1.00002.0000二．雅克比迭代法求解线性方程组Ax = b时，一般当A为低阶稠密矩阵时，用主元消去法求解是常用方法。
但是，对于由工程技术中产生的大型稀疏矩阵方程组(A)阶数很高，但零元素较多，例如求某些偏微分方程数值解所产生的线性方程组)，利用迭代法求解此方程组就是合适的，在计算机内存和运算两方面，迭代法通常都可利用A中有大量零元素的特点。
雅克比迭代法就是众多迭代法中比较早且较简单的一种。
例：用雅克比迭代法求解方程组：10 x- x2=9- x1+10 x2- 2 x3=7-2 x1+10 x2=6编写jacobi函数来实现求解：function x=jacobi(A,b,P,delta,n)%A为n维非奇异阵；b为n维值向量% P为初值；delta为误差界；n为给定的迭代最高次数N=length(b);for k=1:nfor j=1:Nx(j)=(b(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);enderr=abs(norm(x'-P));P=x';if(err<delta)break;endendPx=x';k,err在命令窗口输入：clearA=[10 -1 0;-1 10 -2;0 -2 10];b=[9 7 6]';P=[0 0 0]'; %给出初值x=jacobi(A,b,P,1e-4,20)运行程序，输出如下：P =0.99580.95790.7916k =8err =3.2740e-005x =0.99580.95790.7916三．拉格朗日多项式插值例：已知，，，，090cos 2160cos 2245cos 2330cos ==== 使用Lagrange 多项式插值法计算，，，，，)169cos()78cos()57cos()44cos()35cos( -并给出插值多项式。
编写Lagrange 函数来实现求解：function s=Lagrange(x,y,x0)%Lagrange.m 函数为求已知数据点的拉格朗日插值多项式%x 为数据点的x 坐标向量%y 为数据点的y 坐标向量%x0为插值的x 坐标%s 为求得的拉格朗日插值多项式或在x0处的插值syms p;n=length(x); %读取x 向量维数 s=0;for(k=1:n)la=y(k);for(j=1:k-1)la=la*(p-x(j))/(x(k)-x(j));end;for(j=k+1:n)la=la*(p-x(j))/(x(k)-x(j));end;s=s+la;simplify(s);end%对输入参数个数做判断，如果只有两个参数,直接给出插值多项式%如果三个参数 则给出插值点的插值结果 ,第三个参数可以为向量if(nargin==2)s=subs(s,'p','x');s=collect(s); %多项式展开s=vpa(s,4); %把系数取到6位精度表达elsem=length(x0); %读取x0长度%分别对x0的每一个分量插值for i=1:mtemp(i)=subs(s,'p',x0(i));end%得到的是系列插值点的插值结果s=temp;end在命令窗口输入：clearx=[pi/4, pi/6,pi/3,pi/2];y=[cos(pi/4), cos(pi/6), cos(pi/3), cos(pi/2)];x0=[-35*pi/180, 44*pi/180, 57*pi/180,78*pi/180, 169*pi/180];disp('角度')du=[-35 44 57 78 169]disp('插值结果')yt=Lagrange(x,y,x0)disp('cos 函数值')yreal=[cos(-35*pi/180)cos(44*pi/180)cos(57*pi/180)cos(78*pi/180)cos(169*pi/180)]'disp('插值与函数值误差')dy=yt-yrealyt= Lagrange(x,y)运行程序，输出如下：角度du =-35 44 57 78 169插值结果yt =0.6394 0.7194 0.5446 0.2086 -1.0890cos 函数值yreal =0.8192 0.7193 0.5446 0.2079 -0.9816插值与函数值误差dy =-0.1798 0.0000 -0.0001 0.0006 -0.1074yt =0.1365*x^3 - 0.6731*x^2 + 0.09636*x + 0.9805四．埃尔米特插值法埃尔米特插值要求插值多项式在插值点处的取值与函数值相等，而且还要求导数相同。
例：根据下表所列数据点求出其Hermite插值多项式，并计算当x=2.0时的编写Hermite函数来实现求解：function f = Hermite(x,y,y_1,x0)%Hermite.m插值求已知数据点的埃尔米特插值多项式%x为数据点的x坐标向量%y为数据点的y坐标向量%x0为插值的x坐标% y_1为函数的导数向量%f为求得的埃尔米特插值多项式或在x0处的插值syms t;f = 0.0;if(length(x) == length(y))if(length(y) == length(y_1))n = length(x);elsedisp('y和y的导数的维数不相等！');return;endelsedisp('x和y的维数不相等！');return;endfor i=1:nh = 1.0;a = 0.0;for j=1:nif( j ~= i)h = h*(t-x(j))^2/((x(i)-x(j))^2);a = a + 1/(x(i)-x(j));endendf = f + h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));if(i==n)if(nargin == 4)f = subs(f,'t',x0);elsef = vpa(f,6);endendend在命令窗口输入：clearx=[ 1 1.2 1.4 1.6 1.8];y=[1 1.0954 1.1832 1.2649 1.3416];y_1=[0.5000 0.4564 0.4226 0.3953 0.3727];disp('显示Hermite 插值多项式：')f = Hermite (x,y,y_1)disp('显示在 x=2处的Hermite 插值：')f = Hermite (x,y, y_1,2)运行程序，输出如下：显示Hermite 插值多项式：f =24414.1*(t - 1.8)^2*(t - 1.6)^2*(t - 1.2)^2*(t - 1.0)^2*(0.4226*t + 0.59156) - 678.168*(27.5773*t - 50.9807)*(t - 1.6)^2*(t - 1.4)^2*(t - 1.2)^2*(t - 1.0)^2 - 10850.7*(10.1455*t - 17.4978)*(t - 1.8)^2*(t - 1.4)^2*(t - 1.2)^2*(t - 1.0)^2 + 678.168*(21.3333*t - 20.3333)*(t - 1.8)^2*(t - 1.6)^2*(t - 1.4)^2*(t - 1.2)^2 + 10850.7*(9.58473*t - 10.4063)*(t - 1.8)^2*(t - 1.6)^2*(t - 1.4)^2*(t - 1.0)^2显示在 x=2处的Hermite 插值：f =1.4112五．复化梯形公式 例：用复化梯形积分公式求函数dx x5.15.1-2)-11（的积分，取区间个数n 为10。