数值线性代数第二版徐树方高立张平文上机习题第三章实验报告材料
- 格式:doc
- 大小:130.00 KB
- 文档页数:7
实验报告实验项目名称函数逼近与快速傅里叶变换实验室数学实验室所属课程名称数值逼近实验类型算法设计实验日期班级学号姓名成绩512*x^10 - 1280*x^8 + 1120*x^6 - 400*x^4 + 50*x^2 - 1并得到Figure,图像如下:实验二:编写程序实现[-1,1]上n阶勒让德多项式,并作画(n=0,1,…,10 在一个figure中)。
要求:输入Legendre(-1,1,n),输出如a n x n+a n-1x n-1+…多项式。
在MATLAB的Editor中建立一个M-文件,输入程序代码,实现勒让德多项式的程序代码如下:function Pn=Legendre(n,x)syms x;if n==0Pn=1;else if n==1Pn=x;else Pn=expand((2*n-1)*x*Legendre(n-1)-(n-1)*Legendre(n-2))/(n);endx=[-1:0.1:1];A=sym2poly(Pn);yn=polyval(A,x);plot (x,yn,'-o');hold onend在command Windows中输入命令:Legendre(10),得出的结果为:Legendre(10)ans =(46189*x^10)/256 - (109395*x^8)/256 + (45045*x^6)/128 - (15015*x^4)/128 + (3465*x^2)/256 - 63/256并得到Figure,图像如下:实验三:利用切比雪夫零点做拉格朗日插值,并与以前拉格朗日插值结果比较。
在MATLAB的Editor中建立一个M-文件,输入程序代码,实现拉格朗日插值多项式的程序代码如下:function [C,D]=lagr1(X,Y)n=length(X);D=zeros(n,n);D(:,1)=Y';for j=2:nfor k=j:nD(k,j)=(D(k,j-1)- D(k-1,j-1))/(X(k)-X(k-j+1));endendC=D(n,n);for k=(n-1):-1:1C=conv(C,poly(X(k)));m=length(C);C(m)= C(m)+D(k,k);end在command Windows 中输入如下命令:clear,clf,hold on;k=0:10;X=cos(((21-2*k)*pi)./22); %这是切比雪夫的零点Y=1./(1+25*X.^2);[C,D]=lagr1(X,Y);x=-1:0.01:1;y=polyval(C,x);plot(x,y,X,Y,'.');grid on;xp=-1:0.01:1;z=1./(1+25*xp.^2);plot(xp,z,'r')得到Figure ,图像如下所示:比较后发现,使用切比雪夫零点做拉格朗日插值不会发生龙格现象。
便可求得2.设为两个上三角矩阵,而且线性方程组是非奇异的,试给出一种运算量为的算法,求解该方程组。
因,故为求解线性方程组,可先求得上三角矩阵T的逆矩阵,依照上题的思想我们很容易得到计算的算法。
于是对该问题我们有(1)计算上三角矩阵T的逆矩阵,算法如下:算法 1(求解上三角矩阵的逆矩阵,回代法。
该算法的的运算量为)(2)计算上三角矩阵。
运算量大约为.(3)用回代法求解方程组:.运算量为;(4)用回代法求解方程组:运算量为。
算法总运算量大约为:3.证明:如果是一个Gauss变换,则也是一个Gauss变换。
按Gauss变换矩阵的定义,易知矩阵是Gauss变换。
下面我们只需证明它是Gauss变换的逆矩阵。
事实上注意到,则显然有从而有比较比较向量和可以发现Gauss变换L应具有功能:使向量的第二行加上第一行的2倍;使向量的第三行加上第一行的2倍。
于是Gauss变换如下5.证明:如果有三角分解,并且是非奇异的,那么定理1,其中都是单位下三角阵,都是上三角阵。
因一个单位下三角阵与一个上三角阵相等,故此,它们都必是单位矩阵。
即,6.设的定义如下有满足的三角分解。
[证明] 令 是单位下三角阵,是上三角阵。
定义如下容易验证:7.设A对称且,并假定经过一步Gauss消去之后,证明仍是对称阵。
其中,将A分块为由A的对称性,对称性则是显而易见的。
8.设是严格对角占优阵,即试证:矩阵仍是严格对角占优阵。
由此推断:对于对称的严格对角占优矩阵来说,于是主对角线上的元素满足 (1)非主对角线上的元素满足由于A是严格对角占优的,即故 (2)即,矩阵仍是严格对角占优阵。
9.设有三角分解。
指出当把Gauss消去法应用于矩阵时,怎样才为上三角矩阵U。
而这一组的初等行变换对应的变换矩阵就是,即这就是说,方程组和是同解方程。
而后者是上三角形方程组,可运用本章2求解。
这样我们就不必存储L,通求解方程组,来求解原方程组(1)用初等变换化;(2)利用回代法求解方程组。
一.实验任务用MA TLAB 语言编写连续函数最佳平方逼近的算法程序(函数式M 文件)。
并用此程序进行数值试验,写出实验报告。
二.实验方法最佳平方逼近方法采用基于正交多项式的最佳平方逼近,选择Lengendre 多项式做基。
计算组合系数时,函数的积分采用变步长复化梯形求积法。
三.程序功能和使用说明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. 用M 文件建立数学函数,实现程序通过修改建立数学函数的M 文件以适用不同的被逼近函数。
3.已经考虑一般的情况]1,1[],[)(+-≠∈b a x f ,程序有变量代换的功能。
4.计算组合系数时,函数的积分采用变步长复化梯形求积法5.可根据需要,求出二次、三次、。
最佳平方逼近函数)x s (。
6.最后作出逼近函数)x s (和被逼近函数)(x f 的曲线图可进行比较,分别用绘图函数plot 和fplot 绘图。
7.在matlab 的命令窗口,输入[c,sx]=leastp(@func1,a,b,n),func1是被逼近函数,b 和a 分别是逼近函数的上、下区间,n 为最佳平方逼近的次数,可为任意次数。
四.程序代码(含注释)1. 最佳平方逼近主函数function [c,sx]=leastp(func,a,b,n)%LEASTP.m:least-square fitting with legendre polynomials%func 指被逼近函数,调用需要用句柄%a,b 分别指被逼近函数的区间上下限%n 指最佳平方逼近的次数syms t;syms x;%以Lengendre 多项式为基,构造任意次数的最佳平方逼近多项式p(2)=t;p(1)=1;if n>1for j=3:1:(n+1)p(j)=((2*j-3)*t*p(j-1)-(j-2)*p(j-2))/(j-1);endend%变量代换,区间调整为[-1,1]f=feval(func,(b-a)/2*t+(b+a)/2);%计算组合系数,其中调用变步长复化梯形求积函数trapzfor j=1:1:(n+1)c(j)=(2*j-1)/2*trapz(f*p(j),-1,1);end%将组合系数与对应的最佳平方多项式相乘然后求和,得到最佳逼近函数sx=0;for j=1:1:(n+1)sx=sx+c(j)*p(j);end%将变量替换还原sx=subs(sx,(2*x-a-b)/(b-a));%使用fplot绘制原函数图像f1=feval(func,x);f1=inline(f1);[x,y]=fplot(f1,[a,b]);plot(x,y,'r-','linewidth',1.5);hold on;%使用plot绘制最佳平方逼近函数图像g=linspace(a,b,(b-a)*300);fsx=subs(sx,g);plot(g,fsx,'b-','linewidth',1.5);str=strcat(num2str(n),'次最佳平方逼近');legend('原函数',str);end2. 计算组合系数,变步长复化梯形求积法function To1=trapz(func,a,b)%半分区间复化梯形公式计算定积分%func指需要求积分的原函数%a,b分别指积分上下区间%初值h=b-a;To=(subs(func,a)+subs(func,b))*(b-a)/2;e=1;while e>10^-6%迭代终止条件,前后两次积分值差小于10^-6 H=0;x=a+h/2;while x<bH=H+subs(func,x);%计算出所有二分新出现的值的和x=x+h;endTo1=0.5*(To+h*H);%计算出新的积分值e=abs(To1-To);h=h/2;%继续半分区间,进行迭代计算To=To1;endend3. 以.m文件定义被逼近函数function y=func1(x)y=x*cos(x);end五.实验结果1. 一次最佳平方逼近c =-1.1702 -2.4235sx=1.253290 - 1.211752*x2. 二次最佳平方逼近c =-1.1702 -2.4235 -0.4265sx=-0.159939*x^2 - 0.571997*x + 0.8267873. 三次最佳平方逼近c =-1.1702 -2.4235 -0.4265 1.2216sx=0.381759*x^3 - 2.450495*x^2 + 3.092892*x - 0.3948434. 四次最佳平方逼近c =-1.1702 -2.4235 -0.4265 1.2216 0.3123sx =0.085392*x^4 - 0.301375*x^3 - 0.693864*x^2 + 1.531443*x - 0.082553六.分析与讨论从次数从1到4的最佳平方逼近图像对比可以发现,次数越高,图像拟合效果越好。
数值分析上机实验报告目录1.chapter1舍入误差及有效数 (1)2.chapter2Newton迭代法 (3)3.chapter3线性代数方程组数值解法-列主元Gauss消去法 (7)4.chapter3线性代数方程组数值解法-逐次超松弛迭代法 (8)5.chapter4多项式插值与函数最佳逼近 (10)1.chapter1舍入误差及有效数1.1题目设S N =∑1j 2−1N j=2,其精确值为)11123(21+--N N 。
(1)编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ,计算S N 的通用程序。
(2)编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ,计算S N 的通用程序。
(3)按两种顺序分别计算64210,10,10S S S ,并指出有效位数。
(编制程序时用单精度)(4)通过本次上机题,你明白了什么? 1.2编写相应的matlab 程序 clear;N=input('please input N:'); AValue=((3/2-1/N-1/(N+1))/2); sn1=single(0); sn2=single(0); for i=2:Nsn1=sn1+1/(i*i-1); %从大到小相加的通用程序% endep1=abs(sn1-AValue); for j=N:-1:2sn2=sn2+1/(j*j-1); %从小到大相加的通用程序% endep2=abs(sn2-AValue);fprintf('精确值为:%f\n',AValue);fprintf('从大到小的顺序累加得sn=%f\n',sn1); fprintf('从大到小相加的误差ep1=%f\n',ep1); fprintf('从小到大的顺序累加得sn=%f\n',sn2); fprintf('从小到大相加的误差ep2=%f\n',ep2); disp('================================='); 1.3matlab 运行程序结果 >> chaper1please input N:100 精确值为:0.740050从大到小的顺序累加得sn=0.740049 从大到小相加的误差ep1=0.000001 从小到大的顺序累加得sn=0.740050 从小到大相加的误差ep2=0.000000 >> chaper1please input N:10000 精确值为:0.749900从大到小的顺序累加得sn=0.749852 从大到小相加的误差ep1=0.000048 从小到大的顺序累加得sn=0.749900 从小到大相加的误差ep2=0.000000please input N:1000000精确值为:0.749999从大到小的顺序累加得sn=0.749852 从大到小相加的误差ep1=0.000147 从小到大的顺序累加得sn=0.749999 从小到大相加的误差ep2=0.0000001.4结果分析以及感悟按照从大到小顺序相加的有效位数为:5,4,3。
《数值计算方法》上机实验报告华北电力大学实验名称数值il•算方法》上机实验课程名称数值计算方法专业班级:电力实08学生姓名:李超然学号:200801001008 成绩: 指导教师:郝育黔老师实验日期:2010年04月华北电力大学实验报告数值计算方法上机实验报吿一.各算法的算法原理及计算机程序框图1、牛顿法求解非线性方程*对于非线性方程,若已知根的一个近似值,将在处展开成一阶xxfx ()0, fx ()xkk泰勒公式"f 0 / 2 八八,fxfxfxxxxx 0 0 0 0 0 kkkk2!忽略高次项,有,fxfxfxxx 0 ()()(),,, kkk右端是直线方程,用这个直线方程来近似非线性方程。
将非线性方程的**根代入,即fx ()0, X ,* fxfxxx 0 0 0 0, ,, kkkfx 0 fx 0 0,解出fX 0 *k XX,, k' fx 0 k水将右端取为,则是比更接近于的近似值,即xxxxk, Ik, Ikfx ()k 八XX, Ikk* fx()k这就是牛顿迭代公式。
,2,计算机程序框图:,见,,3,输入变量、输出变量说明:X输入变量:迭代初值,迭代精度,迭代最大次数,\0输出变量:当前迭代次数,当前迭代值xkl,4,具体算例及求解结果:2/16华北电力大学实验报吿开始读入l>k/fx()0?,0fx 0 Oxx,,01* fx ()0XX,,,?10kk, ,1,kN, ?xx, 10输出迭代输出X输出奇异标志1失败标志,3,输入变量、输出变量说明: 结束例:导出计算的牛顿迭代公式,并il •算。
(课本P39例2-16) 115cc (0), 求解结果:10. 75000010.72383710. 72380510. 7238052、列主元素消去法求解线性方程组,1,算法原理:高斯消去法是利用现行方程组初等变换中的一种变换,即用一个不为零的数乘 -个 方程后加只另一个方程,使方程组变成同解的上三角方程组,然后再自下而上 对上三角3/16华北电力大学实验报告方程组求解。
燕山大学课程设计题目:Gauss消去法及其变形方法求解线性方程组学院(系):年级专业:学号:学生姓名:指导教师:教师职称:燕山大学课程设计(论文)任务书院(系):理学院基层教学单位:信息与计算科学系年月日燕山大学课程设计评审意见表目录摘要 (2)一、问题的提出 (2)二、数学模型的建立 (3)三、计算机程序 (3)四、计算结果及分析 (5)五、参考文献 (20)六、心得体会 (20)Gauss 消去法是求解线性方程组的最基本的方法之一。
实际计算中,线性方程组的系数矩阵常常具有对称正定性,即其各阶顺序主子式及全部特征值均大于零。
矩阵的这一特性使它的三角分解也有更简单的形式,从而导出一些特殊的解法,如平方根法和改进的平方根法。
关键词:线性方程组,Gauss 消去法,平方根法,改进的平方根法,Hilbert 矩阵。
一、问题的提出1、先用你所熟悉的计算机语言将不选主元和列主元Gauss 消去法编写成通用程序,然后用你编写的程序求解下面的84阶方程组⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛141515151576816816816816816848382321 x x x x x x最后,将你的计算结果与方程组的精确解进行比较,并就此谈谈你对Gauss 消去法的看法。
2、先用你所熟悉的计算机语言将平方根法和改进的平方根法编写成通用程序,然后用你编写的程序求解对称正定方程组b Ax =,其中 (1)b 随机地选取,系数矩阵为100阶矩阵⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛1011101110111011101110(2)系数矩阵为40阶Hilbert 矩阵,即系数矩阵A 的第i 行第j 列元素为11-+=j i a ij ,向量b 的第i 个分量为∑=-+=nj i j i b 111。
3、用1中的程序求解2中的两个方程组并比较所有的计算结果,然后评论各个方法的优劣。
。
二、数学模型的建立1、Gauss 消去法:用消去法解方程组的基本思想是用逐次消去未知数的方法把原来方程组Ax=b 化为与其等价的三角方程组,而求解三角方程组就容易了。
数值分析上机实验报告实验报告插值法与数值积分实验(数值计算方法,3学时)一实验目的1.掌握不等距节点下的牛顿插值公式以及拉格朗日插值公式。
2.掌握复化的梯形公式、辛扑生公式、牛顿-柯特斯公式计算积分。
3. 会用龙贝格公式和高斯公式计算积分。
二实验内容用拉格朗日插值公式计算01.54.1==y x 以及所对应的近似值。
用牛顿插值公式求)102(y 的近似值。
三实验步骤(算法)与结果1拉格朗日插值法:(C 语言版)#include "Stdio.h" #include "Conio.h"int main(void) {float X[20],Y[20],x; int n;void input(float *,float *,float *,int *); float F(float *,float *,float,int); input(X,Y,&x,&n);printf("F(%f)=%f",x,F(X,Y,x,n));getch(); return 0; }void input(float *X,float *Y,float *x,int *n) {int i;printf("Please input the number of the data:");scanf("%d",n);printf("\nPlease input the locate of each num:\n");for(i=0;i<*n;i++){scanf("%f,%f",X+i,Y+i);}printf("\nPlease input the chazhi:"); scanf("%f",x);}float F(float *X,float *Y,float x,int n){int i,j;float Lx,Fx=0;for(i=0;i<n;i++)< p="">{Lx=1;for(j=0;j<n;j++)< p="">{if(j!=i) Lx=Lx*((x-*(X+j))/(*(X+i)-*(X+j))); } Fx=Fx+Lx*(*(Y+i));}return Fx;}得出结果如图:所以Y(1.4)=3.7295252#include#define N 10double X[N], Y[N], A[N][N];int n;double Newton(double x);double f(double x);void main() {printf("请输入已知x与对应y=f(x)的个数: n = "); scanf("%d", &n);getchar();if(n>N||n<=0) {printf("由于该维数过于犀利, 导致程序退出!"); return;}printf("\n请输入X[%d]: ", n);for (int i=0; i<="" p="">scanf("%lf", &X[i]);getchar();printf("\n请输入Y[%d]: ", n);for (i=0; i<="" p="">scanf("%lf", &Y[i]);getchar();double x;printf("\n请输入所求结点坐标x = ");scanf("%lf", &x);getchar();printf("\nf(%.4lf)≈%lf\n\n", x, Newton(x));}double Newton(double x) {int i, j;// 求均差for (i =0; i<="" p="">A[i][0] = Y[i];for (i=1; i<="" p="">for (j =1; j<=i; j++)A[i][j] = (A[i][j-1] - A[i-1][j-1]) / (X[i] - X[i-j]); // 求结点double result = A[0][0];for (i=1; i<="">double tmp = 1.0;for (int j=0; j<="" p="">tmp *= (x - X[j]);result += tmp * A[i][i];}return result;}四实验收获与教师评语</n;j++)<></n;i++)<>。
数值分析上机实践报告班级:计算机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<<"次。
.. .下载可编辑. 第三章上机习题
用你所熟悉的的计算机语言编制利用QR分解求解线性方程组和线性最小二乘问题的通
用子程序,并用你编制的子程序完成下面的计算任务: (1)求解第一章上机习题中的三个线性方程组,并将所得的计算结果与前面的结果相比较,说明各方法的优劣;
(2)求一个二次多项式+bt+cy=at2,使得在残向量的2数下最小的意义下拟合表3.2中的数据; 表 3.2 ti -1 -0.75 -0.5 0 0.25 0.5 0.75
yi 1 0.8125 0.75 1 1.3125 1.75 2.3125
(3)在房产估价的线性模型 111122110xaxaxaxy
中,1121,,,aaa分别表示税、浴室数目、占地面积、车库数目、房屋数目、居室数目、房龄、建筑类型、户型及壁炉数目,y代表房屋价格。现根据表3.3和表3.4给出的28组数据,求出模型中参数的最小二乘结果。 (表3.3和表3.4见课本P99-100) 解 分析:
(1)计算一个Householder变换H:
由于TTvvIwwIH2,则计算一个Householder变换H等价于计算相应的v、。
其中)/(2,||||12vvexxvT。 在实际计算中,
为避免出现两个相近的数出现的情形,当01x时,令212221||||)(-xxxxvn;
为便于储存,将v规格化为1/vvv,相应的,变为)/(221vvvT 为防止溢出现象,用||||/xx代替 (2)QR分解: 利用Householder变换逐步将nmAnm,转化为上三角矩阵AHHHnn11,则有 .. .下载可编辑.
0RQA,其中nHHHQ21,:),:1(nR。
在实际计算中,从nj:1,若mj,依次计算)),:((jmjAx对应的)1()1()~(kmkmjH 即对应的jv,j,将)1:2(jmvj储存到),:1(jmjA,j储存到)(jd,迭代结束
后再次计算Q,有~001jjjHIH,nHHHQ21(mn时1-21nHHHQ) (3)求解线性方程组bAx或最小二乘问题的步骤为 i 计算A的QR分解;
ii 计算bQcT11,其中):1(:,1nQQ
iii 利用回代法求解上三角方程组1cRx (4)对第一章第一个线性方程组,由于R的结果最后一行为零,故使用前代法时不计最后一行,而用运行结果计算84x。
运算matlab程序为 1 计算Householder变换 [v,belta]=house(x) function [v,belta]=house(x) n=length(x); x=x/norm(x,inf); sigma=x(2:n)'*x(2:n); v=zeros(n,1); v(2:n,1)=x(2:n); if sigma==0 belta=0; else alpha=sqrt(x(1)^2+sigma); if x(1)<=0 v(1)=x(1)-alpha; else v(1)=-sigma/(x(1)+alpha); end belta=2*v(1)^2/(sigma+v(1)^2); v=v/v(1,1); end end .. .下载可编辑. 2 计算A的QR分解 [Q,R]=QRfenjie(A) function [Q,R]=QRfenjie(A) [m,n]=size(A); Q=eye(m); for j=1:n if j [v,belta]=house(A(j:m,j)); H=eye(m-j+1)-belta*v*v'; A(j:m,j:n)=H*A(j:m,j:n); d(j)=belta; A(j+1:m,j)=v(2:m-j+1); end end R=triu(A(1:n,:)); for j=1:n if j H=eye(m); temp=[1;A(j+1:m,j)]; H(j:m,j:m)=H(j:m,j:m)-d(j)*temp*temp'; Q=Q*H; end end end
3 解下三角形方程组的前代法 x=qiandaifa(L,b) function x=qiandaifa(L,b) n=length(b); for j=1:n-1 b(j)=b(j)/L(j,j); b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j); end b(n)=b(n)/L(n,n); x=b; end
4 求解第一章上机习题中的三个线性方程组 ex3_1
clear;clc; %第一题 A=6*eye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1); b=[7;15*ones(82,1);14]; n=length(A); .. .下载可编辑. %QR分解 [Q,R]=QRfenjie(A); c=Q'*b; x1=huidaifa(R(1:n-1,1:n-1),c(1:n-1)); x1(n)=c(n)-R(n,1:n-1)*x1; %不选主元Gauss消去法 [L,U]=GaussLA(A); x1_1=Gauss(A,b,L,U); %列主元Gauss消去法 [L,U,P]=GaussCol(A); x1_2=Gauss(A,b,L,U,P); %解的比较 figure(1); subplot(1,3,1);plot(1:n,x1);title('QR分解'); subplot(1,3,2);plot(1:84,x1_1);title('Gauss'); subplot(1,3,3);plot(1:84,x1_2);title('PGauss'); %第二题第一问 A=10*eye(100)+diag(ones(1,99),-1)+diag(ones(1,99),1); b=round(100*rand(100,1)); n=length(A); %QR分解 tic;[Q,R]=QRfenjie(A); c=Q'*b; x2=huidaifa(R,c);toc; %不选主元Gauss消去法 tic;[L,U]=GaussLA(A); x2_1=Gauss(A,b,L,U);toc; %列主元Gauss消去法 tic;[L,U,P]=GaussCol(A); x2_2=Gauss(A,b,L,U,P);toc; %平方根法 tic;L=Cholesky(A); x2_3=Gauss(A,b,L,L');toc; %改进的平方根法 tic;[L,D]=LDLt(A); x2_4=Gauss(A,b,L,D*L');toc; %解的比较 figure(2); subplot(1,5,1);plot(1:n,x2);title('QR分解'); subplot(1,5,2);plot(1:n,x2_1);title('Gauss'); subplot(1,5,3);plot(1:n,x2_2);title('PGauss'); subplot(1,5,4);plot(1:n,x2_3);title('平方根法'); subplot(1,5,5);plot(1:n,x2_4);title('改进的平方根法'); %第二题第二问 .. .下载可编辑. A=hilb(40); b=sum(A); b=b'; n=length(A); [Q,R]=QRfenjie(A); c=Q'*b; x3=huidaifa(R,c); %不选主元Gauss消去法 [L,U]=GaussLA(A); x3_1=Gauss(A,b,L,U); %列主元Gauss消去法 [L,U,P]=GaussCol(A); x3_2=Gauss(A,b,L,U,P); %平方根法 L=Cholesky(A); x3_3=Gauss(A,b,L,L'); %改进的平方根法 [L,D]=LDLt(A); x3_4=Gauss(A,b,L,D*L'); %解的比较 figure(3); subplot(1,5,1);plot(1:n,x3);title('QR分解'); subplot(1,5,2);plot(1:n,x3_1);title('Gauss'); subplot(1,5,3);plot(1:n,x3_2);title('PGauss'); subplot(1,5,4);plot(1:n,x3_3);title('平方根法'); subplot(1,5,5);plot(1:n,x3_4);title('改进的平方根法');
5 求解二次多项式 ex3_2
clear;clc; t=[-1 -0.75 -0.5 0 0.25 0.5 0.75]; y=[1 0.8125 0.75 1 1.3125 1.75 2.3125]; A=ones(7,3); A(:,1)=t'.^2; A(:,2)=t'; [Q,R]=QRfenjie(A); Q1=Q(:,1:3); c=Q1'*y'; x=huidaifa(R,c)
6 求解房产估价的线性模型 ex3_3
clear;clc; A=xlsread('E:\temporary\专业课\数值代数\cha3_3_4.xls','A2:L29'); y=xlsread('E:\temporary\专业课\数值代数\cha3_3_4.xls','M2:M29'); [Q,R]=QRfenjie(A);