数值分析MATLAB实验程序.docx

数值分析上机实验报告院系:XXX学号:XXX姓名:XXXX目录一.Gass-Jordan 消去法 1二．雅克比迭代法 2三．拉格朗日多项式插值 4四．埃尔米特插值法 5五．复化梯形公式 7六．龙贝格求积公式 8七．欧拉法 9八．隐式欧拉法 10．Gass-Jordan 消去法Gass-Jordan消去法求解线性方程组的基本思想是将系数矩阵化为对角矩阵, 这样可以直接得到方程组的解x.x2・・x，因此也称为无回代消去法12n例：用列主元Gass-Jordan消去法求解方程组：{2X]+1x2+x3=55X]-1x2+1x3=81X]-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))>max1 max1=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-X]+10 x2- 2 x3=7 -2 x]+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:N x(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

■■ ,cos45 ' ,cos60 ,cos90 0,使用 Lagrange 多2 2 2项式插值法计算cos(-35),cos(44o),osQ7),os(7&),os(169)，并给出插值多项 式编写 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-yreal运行程序，输出如下yt= Lagrange(x,y)-35 445778 169插值结果yt =0.63940.71940.54460.2086-1.0890cos 函数值yreal =0.81920.71930.54460.2079-0.9816角度du =插值与函数值误差dy =-0.17980.0000 -0.00010.0006 -0.1074yt =0.1365*xT - 0.6731*xT + 0.09636*x + 0.9805四．埃尔米特插值法埃尔米特插值要求插值多项式在插值点处的取值与函数值相等，而且还要求导数相同。

例：根据下表所列数据点求出其 Hermite 插值多项式，并计算当 x=2.0 时的 y 值x11.21.41.61.8y11.09541.18321.26491.3416y i0.50000.45640.42260.39530.3727编写 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;end elsedisp('x和y的维数不相等！’)；return； end for i=1:nh = 1.0；a = 0.0； for j=1:n if( j ~= i)h = h*(t-x(j))A2/((x(i)-x(j))A2); 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)；endend end在命令窗口输入：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)A2*(t - 1.4)A2*(t - 1.2)A2*(t - 1.0)人2 - 10850.7*(10.1455*t - 17.4978)*(t - 1.8)人2*亿-1.4)人2*亿-1.2)人2*亿-1.0)人2 + 678.168*(21.3333*t - 20.3333)*(t - 1.8)A2*(t - 1.6)A2*(t - 1.4)A2*(t - 1.2)A2 + 10850.7*(9.58473*t - 10.4063)*(t - 1.8)A2*(t - 1.6)A2*(t - 1.4)A2*(t - 1.0。