1、数值 计算方法 课程设计 姓名 高振翔 学号 201211011066 班级 信计 12-2 成绩： 1 编写秦九韶算法程序，并用该程序计算多项式 623)( 35 xxxxf 在1.3 1.2, ,1.1x 的值 function value=qinjiushao(A,x) n=length(A); F=zeros(n); F(1)=A(1); for i=1:n-1 F(i+1)=F(i)*x+A(i+1); end value=F(n) disp(真值 ) polyval(A,x) 程序 ： function s=qinjiuzhao(a,x) n=length(a); s=a(1);
2、for k=2:n; s=s*x+a(k); end 3. 调试： a=1 0 3 2 6; s=qinjiuzhao(a,1,1) s=qinjiuzhao(a,1,2) s=qinjiuzhao(a,1,3) 结果： s=9.4035 s=11.2723 s=13.7039 2. *用选列主元高斯消去法 解线性方程组 02 02 0 21 34343232121xxxxxxxxxx计算的 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
3、)*b; %matlab 的计算结果 y= -0.2667 -0.2000 -0.1333 -0.0667 n=length(b); x=zeros(n,1); %方程个数 n;x 未知向量 % 以下 消 为 去过程 for k=1:n-1 % if A(k,k)=0; % error(Error); % end for i=k+1:n % A(i,k)=A(i,k)/A(k,k); Aik=A(i,k)/A(k,k) for j=k:n A(i,j)=A(i,j)-Aik*A(k,j); end A b(i)=b(i)-Aik*b(k) end end % 回代 Aik = 0.3333 A
4、 = -3.0000 -1.0000 0 0 0 2.3333 -1.0000 0 0 -1.0000 2.0000 -1.0000 0 0 -1.0000 2.0000 b = 1.0000 -0.3333 0 0 Aik = 0 A = -3.0000 -1.0000 0 0 0 2.3333 -1.0000 0 0 -1.0000 2.0000 -1.0000 0 0 -1.0000 2.0000 b = 1.0000 -0.3333 0 0 Aik = 0 A = -3.0000 -1.0000 0 0 0 2.3333 -1.0000 0 0 -1.0000 2.0000 -1.00
5、00 0 0 -1.0000 2.0000 b = 1.0000 -0.3333 0 0 Aik = -0.4286 A = -3.0000 -1.0000 0 0 0 2.3333 -1.0000 0 0 0 1.5714 -1.0000 0 0 -1.0000 2.0000 b = 1.0000 -0.3333 -0.1429 0 Aik = 0 A = -3.0000 -1.0000 0 0 0 2.3333 -1.0000 0 0 0 1.5714 -1.0000 0 0 -1.0000 2.0000 b = 1.0000 -0.3333 -0.1429 0 Aik = -0.6364
6、 A = -3.0000 -1.0000 0 0 0 2.3333 -1.0000 0 0 0 1.5714 -1.0000 0 0 0 1.3636 b = 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:1 S=b(k); for j=k+1:n S=S-A(k,j)*x(j); end x(k)=S/A(k,k) end x = 0 0 -0.1333 -0.0667 x = 0 -0.2000 -0.1333 -0.0667 x = -0.2667 -0.2000 -0.133
7、3 -0.0667 x x = -0.2667 -0.2000 -0.1333 -0.0667 error=abs(x-ones(n,1) )%误差 error = 1.2667 1.2000 1.1333 1.0667 toc%运行时间命令 运行时间： 348.6710 结构分析： 在用高斯消去法求解方程组的解，化为阶梯型时，主元过小可能产生麻烦，会产生很大的误差， 既小 主元要在分母上，产生的误差变化很大， 所以应避免采用绝对 值最小的主元素，对于一般矩阵来说，最好每一步选取系数矩阵或 消元后的低阶矩阵中绝对值最大的元素作为主元素，使高斯消去法具有较好的稳定性，主要使用列主元消去法！ 小结
8、： 在求解方程组时，使用列主元消去法，先判定方程组的系数矩阵非奇异，然后进行行变换，按列主元消去法化为阶梯型，当计算到系数行列式为 0 时计算停止，然后在回代求解最终求得原方程组的解。 5 *用二分法和 Newton 迭代法求下列方程的正根： ,05.01)1ln ( 22 xxxxx function test clear clc %实验方程： 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
9、,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
10、=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 %- end 9*取 2.0h ，由改进 Euler 格式解微 分方程： 2 22 3 ( 0 ) 5, 0 10y x y y x 计算的 matlab 程序如下： clc,clear tic syms x y f=(2*x*
11、x+3*y*y)1/2); a=0;b=10; h=0.2; n=(b-a)/h+1;%n=(b-a)/h x=0;y=5; szj=x,y; for i=1:n-1%i=1:n y=y+h*subs(f,x,y,x,y); x=x+h; szj=szj;x,y; end szj szj = 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
12、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
13、 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.4401 plot(szj(:,1),szj(:,2),or-) toc e
14、lapsed_time = %cpu 运行时间 194.4530 结构分析：在求解初值问题的常微分方程式有改进的欧拉公式进行求解时具体步骤为：分割求解区间 :已知步长 h=x(k+1)-x(k)=(b-a)/n=0.2， 所以 n=50,区间为（ 0， 10），所以等距剖分a=x0 a=0;b=10;h=0.2; n=(b-a)/h+1; x=0;y=6; szj=x y; for i=1:n-1 l1=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,y; end Warning: Unmatched “end“. (Type “warning off MATLAB:m_warning_end_without_block“ to suppress this warning.) plot(szj(:,1),szj(:,2),dg-) toc elapsed_time = 767.5780 %cpu 运行时间 结构分析： 当求积节点越多时 ，精度越高，显然当节点个数为 1 时即为欧拉法
