科学计算实践精要.docx

题目一：求积分方程的离散解并比较它与真解之间的误差用的复化Sinpson公式求积分方程的近似解，并在区间的一些点上与真解进行比较分析：1.利用复化Sinpson公式替换积分，得到关于x和y的二元方程；2.分别令得到2n+1个关于的线性方程组；3.使用Gauss消去法求解线性方程组,得到这2n+1个点的函数值；4.使用三次样条函数得到数值解；5.画图：在同一图上标出，各处的数值解和解析解，及数值解(利用三次样条构造)和解析解的函数图像。源程序及运行结果：%科学计算实践第1题老杨clear;clc;syms xtn=input(请输入区间分割份数 );a=0; b=1;f=(4*x^3+5*x^2-2*x+5)/(8*(x+1)^2); g=1/(1+t);h=(b-a)/n;xi=a:h/2:b; %将区间分为2n+1个小区间%题设初始条件I=eye(2*n+1);%单位矩阵K=zeros(2*n+1,2*n+1);%K矩阵for i=1:2*n+1for j=1:2*n+1K(i,j)=-xi(i)+subs(g,t,xi(j));endendd=[1,repmat([4,2],1,n-1),4,1]; %Simpson插值的各项系数D=h/6*diag(d);%D权重矩阵B=zeros(2*n+1,1);for i=1:2*n+1 B(i)=subs(f,x,xi(i));endY=linsolve((I-K*D),B);Y1=csapi(xi,Y);%三次样条插值fnplt(Y1);%显示插值结果hold onplot(xi,Y,*r)%插值点%%%%%%%%%%%%%真值检验y2=1/(x+1)^2;%检验函数fplot(1/(x+1)^2,[0,1])%检验函数图像for i=1:2*n+1 Y2(i)=subs(y2,x,xi(i));endplot(xi,Y2,ok)%插值点legend(数值解,数值解插值点,解析解,解析解插值点);n=2时：n=4时：n=40时：Published with MATLAB? R2013a题目二：分别用共轭梯度算法和预处理共轭梯度算法求解线性方程组其中分析：一、共轭梯度法因为矩阵形式为对成正定且对角占优，因此可以使用共轭梯度法求解线性方程组，需要选定初始值，下降方向和下降步长。之前在数值分析课程已经学过共轭梯度法，故对该方法不再过多介绍，详细理论见《数值分析》梅立泉 P92-99页。算法如下：第1步：取初始向量, 计算第k+1步：计算二、预处理共轭梯度法先对矩阵进行适当的预处理，然后再进行共轭梯度法。根据矩阵的实际情况，此处采用对角预优矩阵法。取预优矩阵，M为矩阵A的对角元矩阵。进行预处理后进行计算。源程序及运行结果：主程序：%科学计算实践第2题 老杨%载入矩阵d=zeros(1,20);for i=1:20 d(i)=2^i;endA=diag(d)+diag(-1*ones(19,1),1)+diag(-1*ones(19,1),-1);D=d-2;D(1)=D(1)+1;D(20)=D(20)+1;tic;disp(1. 共轭梯度法解方程);x1=Gongetidu(A,D)toc;tic;disp(2. 调用系统的预处理共轭梯度算法解方程);x2=pcg(A,D)toctic;M=diag(d);disp(3. 调用自编的预处理共轭梯度算法解方程);y=Gongetidu(M^(-1/2)*A*M^(-1/2),M^(-1/2)*D);%预处理x3=linsolve(M^(1/2),y)toc调用的自编共轭梯度算法程序：%共轭梯度法解线性方程组%《数值分析》梅立泉 P99页function x=Gongetidu(A,b,x0,epsa)n=500;if(nargin==2) %如果只输入了A和b，则设置默认的初始搜索方向及误差上限 x0=zeros(size(A,1),1); epsa=1.0e-8;endx=x0;r=b-A*x;d=r;for k=0:(n-1) alpha=(r*r)/(d*A*d); x=x+alpha*d; r2=b-A*x;if ((norm(r2)=epsa)|(k==n-1)) x; k %迭代次数break;end beta=norm(r2)^2/norm(r)^2;d=r2+beta*d;r=r2;end精度为1.0e-6时：1. 共轭梯度法解方程k = 62x1 = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000