西安交大计算方法b大作业

m(k)=(yy(k)-r(k)*m(k+1))/u(k); end %%% 绘制曲线 k=1; nn=100; xx=linspace(0,20,nn); l=0; for for if k=i; break ; else k=i+1; end end h=1; xbar=x(k)-xx(j); xmao=xx(j)-x(k-1); s(j)=(m(k-1)*xbar^3/6+m(k)*xmao^3/6+(y(k-1)-m(k-1)*h^2/6)*xbar+(y(k)m(k)*h^2/6)*xmao)/h; j=1:nn i=2:20 xx(j)<=x(i)
7
');
G
for j=1:n for i=1:m G(i,j)=x(1,i)^(j-1); end end %%% 建立矩阵 for k=1:n for i=k:m c(k)=G(i,k)^2+c(k); end c(k)=-sign(G(k,k))*(c(k)^0.5); w(k)=G(k,k)-c(k); % 建立 w 来 存放 ω for j=k+1:m w(j)=G(j,k); end b(k)=c(k)*w(k); %%% 变换矩阵 G(k,k)=c(k); for j=k+1:n+1 q=0; for i=k:m q=w(i)*G(i,j)+q; end s=q/b(k); for i=k:m G(i,j)=s*w(i)+G(i,j); end end end %%% 求解三角方程 a(n)=G(n,n+1)/G(n,n); for i=n-1:(-1):1 for j=i+1:n f=G(i,j)*a(j)+f; end a(i)=(G(i,n+1)-f)/G(i,i); f=0; end a %%% 回代过程 p=zeros(1,m); for j=1:m for i=1:n p(j)=p(j)+a(i)*x(j)^(i-1); end end %a（ i ）存放各级系数 Rx=h1 Gk-1 到 G k Qk
5
sp(j)=-m(k-1)*(x(k)-xx(j))^2/(2*h)+m(k)*(xx(j)-x(k-1))^2/(2*h)+(y(k)y(k-1))/h-(m(k)-m(k-1))*h/6; l(j+1)=(1+sp(j)^2)^0.5*(20/nn)+l(j); %求解光缆长度 end %%% 绘图 plot(xx,s, disp([ 'r-' , 'linewidth' ,1.5) '? ×＇ ])
8
plot(x,p, 'r*' ,x,p, '-' ); E2=0; % 用 E2 来存放误差 %%% 误差求解 for i=n+1:m E2=G(i,n+1)^2+E2; end E2=E2^0.5; disp( '误差 为 '); disp(E2); t=0; for i=1:m t=t+p(i); end t=t/m; %%% 平均温度 disp([ '平均温度为 ',num2str(t), '℃ '])
实现思想： 此题中所给数据点数目较多， 采用拉格朗日插值法或者牛顿插值法需 要很高次的多项式，计算困难，误差大；采用样条插值计算量虽然不大，但是存 放参数 Mi 的量很大，且没有一个统一的数学公式来表示，也不是很方便。所以 可考虑用最小二乘法进行拟合。 计算过程中， 分别使用二次函数、 三次函数以及 四次函数，计算其相应的系数，估算误差并作图比较各个函数之间的区别。 算法结构：（ 参考课本 P123） 1.1[ 形成矩阵 Qk] 1.2[ 变换 Gk-1 到 Gk] 2.[ 求解三角方程 ] 3.[ 计算误差 ] 源代码：
实验结果： 二次函数拟合，结果如下图所示
图 3-1 二次函数拟合结果
9
三次函数拟合，结果如下图所示
图 3-2 三次函数拟合结果 四次函数拟合，结果如下图所示
图 3-3 四次函数拟合结果
10
结果对比 ： 将二次函数、三次函数和四次函数拟合结果绘制在同一个坐标内， 如图 3-4 所示。其计算误差结果见表 3-input( e=10^(-4); Nmax=1000; for ' 输入初值 ' ); %误差设定 %迭代最大次数限定
n=1:Nmax
f0=f2(x1); if abs(f2(x1))<e ' 输出的 f(x) 已经足够小 ' );
fprintf( x=x1; break else F0=f3(x1); x=x1-f0/F0; if
s0=4/(16^n*(8*n+1))-2/(16^n*(8*n+4))-1/(16^n*(8*n+5))-1/(16^n*(8*n+6) ); n=n+1; end %分别求解各项并求和 for k=n-1:-1:0 a1=4/(16^k*(8*k+1)); a2=2/(16^k*(8*k+4)); a3=1/(16^k*(8*k+5)); a4=1/(16^k*(8*k+6)); s1=a1+s1; s2=a4+a3+a2+s2; end S=vpa(s1-s2,m)
《计算方法 B》上机 实验报告
学院： 班级 : 姓名： 学号：
机械工程学院
2015 年 12 月 22 日
1
1. 计算以下和式： S
n 0
1 16
n
4 8n 1 8n
2 4
1 8n 5
1 8n 6
，要求：
(1) 若保留 11 个有效数字，给出计算结果，并评价计算的算法； (2) 若要保留 30 个有效数字，则又将如何进行计算。 实现思想：
20 19
P110；
20
l =
0
1 + (f( x )) dx
k 0
'
2
k 1 k
1 + (f( x )) dx
0
'
2
源程序：
clear; clc; x=0:20; y=[9.01 8.96 7.96 7.97 8.02 9.05 10.13 11.18 12.26 13.28 13.32 12.61 11.29
4
d(i)=6*d(i+1); end d(1)=0; d(21)=0; %%% 带状矩阵求解（追赶法） a=0.5*ones(1,21); b=2*ones(1,21); c=0.5*ones(1,21); a(1)=0; c(21)=0; u=ones(1,21); u(1)=b(1); r=c; yy(1)=d(1); %%% 追 for k=2:21 l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*r(k-1); yy(k)=d(k)-l(k)*yy(k-1); end %%% 赶 m(21)=yy(21)/u(21); for k=20:-1:1
5
2
) .
, 带入不同的初值即可解出不同的根
源代码 :
function y=f2(x) %定义原函数 y=6*x.^5-45*x.^2+20; function y=f3(x) % 定义原函数倒数
y=30*x^4-90*x;
i=-5:0.1:5; y=f2(i); plot(i,y) hold plot(i,0, on '-' ) %画出原函数图像
图 3-4 拟合结果对比分析
11
4. 设计算法，求出非线性方程 6 x5 45 x2 20 0 的所有实根，并使 误差不超过 10 4 。
算法思想： 本题可采用牛顿法迭代求解 , 令 f(x) x k+1 = x k 根据函数图像可以找出根的大致分布区间 f(x k ) f(x
' k
= 6 x - 45 x + 20 , 得带格式为
2
实验结果： 11 位有效数字计算结果如图 1 所示； 30 为有效数字计算结果如图 2 所 示。
图 1.11 位有效数字计算结果
图 2.30 为有效数字计算结果
3
1. 某通信公司在一次施工中， 需要在水面宽度为 20 米的河沟底部沿 直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形 进行初步探测， 从而估计所需光缆的长度， 为工程预算提供依据。 已探测到一组等分点位置的深度数据
以上问题出现了近似数相减的问题， 为了减小误差， 可分别求得减数之和以 及被减数之和，最后将两者相减。另外，减数与被减数求和均为同号计算，按照 绝对值递增顺序相加可减小舍入误差。 此题中对有效数字有要求， 因而计算时首
先需要根据有效数字位数计算得出迭代次数，以保证计算值的精度。
源程序：
m=input( ' 输入有效数字个数 m=' ); s0=1;s1=0;s2=0;n=0; %判断迭代次数 while s0>=0.5*10^-(m-1)
clear;clc; x=0:24; y=[15 14 14 14 14 15 16 18 20 20 23 25 28 31 34 31 29 27 25 24 22 20 18 17 16]; m=length(x); n=input( ' 请输入函数的次数 plot(x,y, 'k.' ,x,y,'-' ) grid; hold on ; n=n+1; G=zeros(m,n+1); G(:,n+1)=y'; c=zeros(1,n); %建立 c 来存 放 σ q=0; f=0; b=zeros(1,m); % 建 立 b 用来存放 β%%% 形成矩阵