数值上机作业.doc

数值计算方法课程设计姓名 高振翔 学号201211011066班级 信计12-2 成绩：1编写秦九韶算法程序，并用该程序计算多项式在的值function value=qinjiushao(A,x) n=length(A); F=zeros(n); F(1)=A(1); for i=1:n-1F(i+1)=F(i)*x+A(i+1); endvalue=F(n) disp(真值) polyval(A,x)程序：functions=qinjiuzhao(a,x) n=length(a); s=a(1); fork=2:n; s=s*x+a(k); end 3. 调试： a=10326; s=qinjiuzhao(a,1,1) s=qinjiuzhao(a,1,2) s=qinjiuzhao(a,1,3) 结果：s=9.4035s=11.2723 s=13.70392. *用选列主元高斯消去法解线性方程组 计算的matlab程序：tic%运行时间命令A=-3 -1 0 0;-1 2 -1 0;0 -1 2 -1;0 0 -1 2; b=1 0 0 0; ;%A系数矩阵，b为n维向量y=inv(A)*b; %matlab的计算结果y= -0.2667 -0.2000 -0.1333 -0.0667n=length(b); x=zeros(n,1); %方程个数n;x未知向量%以下消为去过程for k=1:n-1% if A(k,k)=0;% error(Error);% endfor i=k+1:n% A(i,k)=A(i,k)/A(k,k); Aik=A(i,k)/A(k,k)for j=k:nA(i,j)=A(i,j)-Aik*A(k,j);endAb(i)=b(i)-Aik*b(k)endend%回代Aik =0.3333A = -3.0000 -1.0000 0 0 0 2.3333 -1.0000 0 0 -1.0000 2.0000 -1.0000 0 0 -1.0000 2.0000b = 1.0000 -0.3333 0 0Aik =0A = -3.0000 -1.0000 0 0 0 2.3333 -1.0000 0 0 -1.0000 2.0000 -1.0000 0 0 -1.0000 2.0000b = 1.0000 -0.3333 0 0Aik = 0A = -3.0000 -1.0000 0 0 0 2.3333 -1.0000 0 0 -1.0000 2.0000 -1.0000 0 0 -1.0000 2.0000b = 1.0000 -0.3333 0 0Aik = -0.4286A = -3.0000 -1.0000 0 0 0 2.3333 -1.0000 0 0 0 1.5714 -1.0000 0 0 -1.0000 2.0000b = 1.0000 -0.3333 -0.1429 0Aik = 0A = -3.0000 -1.0000 0 0 0 2.3333 -1.0000 0 0 0 1.5714 -1.0000 0 0 -1.0000 2.0000b = 1.0000 -0.3333 -0.1429 0Aik = -0.6364A = -3.0000 -1.0000 0 0 0 2.3333 -1.0000 0 0 0 1.5714 -1.0000 0 0 0 1.3636b = 1.0000 -0.3333 -0.1429 -0.0909 x(n)=b(n)/A(n,n)x = 0 0 0 -0.0667 for k=n-1:-1:1S=b(k); for j=k+1:nS=S-A(k,j)*x(j);endx(k)=S/A(k,k)endx = 0 0 -0.1333 -0.0667x = 0 -0.2000 -0.1333 -0.0667x = -0.2667 -0.2000 -0.1333 -0.0667 xx = -0.2667 -0.2000 -0.1333 -0.0667 error=abs(x-ones(n,1) )%误差error = 1.2667 1.2000 1.13331.0667toc%运行时间命令运行时间： 348.6710结构分析：在用高斯消去法求解方程组的解，化为阶梯型时，主元过小可能产生麻烦，会产生很大的误差，既小主元要在分母上，产生的误差变化很大，所以应避免采用绝对值最小的主元素，对于一般矩阵来说，最好每一步选取系数矩阵或消元后的低阶矩阵中绝对值最大的元素作为主元素，使高斯消去法具有较好的稳定性，主要使用列主元消去法！小结：在求解方程组时，使用列主元消去法，先判定方程组的系数矩阵非奇异，然后进行行变换，按列主元消去法化为阶梯型，当计算到系数行列式为0时计算停止，然后在回代求解最终求得原方程组的解。5*用二分法和Newton迭代法求下列方程的正根： function testclearclc%实验方程：3*x.2+x+2*exp(x)=0%原函数f=(x)3*x.2+x-2*exp(x);%导函数df=(x)6*x+1-2*exp(x);%原函数-1 0上图像（有根范围）fplot(f,-1 0)hold on%牛顿切线法x1,n1=fnewton(f,df,-0.5);disp(sprintf(牛顿切线法n%f附近根：%fn迭代次数：%d,-0.5,x1,n1)%二分法x2,n2=f2fen(f,-1,0);disp(sprintf(二分法n%f,%f上根：%fn迭代次数：%d,-1,0,x2,n2)plot(x1,f(x1),xr,x2,f(x2),+g)%-牛顿切线法- function x,n=fnewton(f,df,x0) x=x0;%初值 delta=1; n=0;%迭代次数下同 while abs(delta)1e-6 delta=f(x)/df(x); x=x-delta; n=n+1; end end%-二分法- function x,n=f2fen(f,a,b) xab=a;b;%两端点值 pab=sign(f(xab); n=0; while diff(xab)1e-6 x=mean(xab); p=sign(f(x); n=n+1; if p,break;end xab(p=pab)=x; end end%-end9*取，由改进Euler格式解微分方程： 计算的matlab程序如下：clc,clearticsyms x yf=(2*x*x+3*y*y)1/2);a=0;b=10;h=0.2;n=(b-a)/h+1;%n=(b-a)/hx=0;y=5;szj=x,y;for i=1:n-1%i=1:ny=y+h*subs(f,x,y,x,y);x=x+h;szj=szj;x,y;endszjszj =1.0e+007 * 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0000 0.0001 0.0000 0.0001 0.0000 0.0001 0.0000 0.0002 0.0000 0.0003 0.0000 0.0003 0.0000 0.0005 0.0000 0.0006 0.0000 0.0008 0.0000 0.0011 0.0000 0.0015 0.0000 0.0021 0.0000 0.0028 0.0000 0.0038 0.0000 0.0051 0.0000 0.0068 0.0000 0.0092 0.0000 0.0123 0.0000 0.0166 0.0000 0.0224 0.0000 0.0301 0.0000 0.0406 0.0000 0.0546 0.0000 0.0736 0.0000 0.0990 0.0000 0.1333 0.0000 0.1795 0.0000 0.2417 0.0000 0.3255 0.0000 0.4382 0.0000 0.5900 0.0000 0.7944 0.0000 1.0696 0.0000 1.4401plot(szj(:,1),szj(:,2),or-)tocelapsed_time = %cpu运行时间194.4530结构分析：在求解初值问题的常微分方程式有改进的欧拉公式进行求解时具体步骤为：分割求解区间:已知步长h=x(k+1)-x(k)=(b-a)/n=0.2，所以n=50,区间为（0，10），所以等距剖分a=x0x1x2 a=0;b=10;h=0.2; n=(b-a)/h+1; x=0;y=6; szj=x y; for i=1:n-1l1=subs(f,x,y,x,y); l2=subs(f,x,y,x+h/2,y+l1*h/2); l3=subs(f,x,y,x+h/2,y+l2*h/2); l4=subs(f,x,y,x+h,y+l3*h); y=y+h*(l1+2*l2+2*l3+l4)/6; x=x+h; szj=szj;x