北科大 计算方法matlab作业

北京科技大学 冶金与生态工程学院 冶金工程

冶研1201班s20120273

2012级研究生 《计算方法》作业

姓名：

学号： s20120273 专业： 冶金工程 学院： 冶金与生态工程学院 成绩：__________________ 任课教师：数理学院 丁军

2012年11月20日

18811349978 1

北京科技大学 冶金与生态工程学院 冶金工程

实验一 牛顿下山法

实验目的

1. 掌握牛顿下山法求解方程根的推导原理。 2. 理解牛顿下山法的具体算法与相应程序的编写。 实验内容

采用牛顿下山法求方程2x3-5x-17=0在2附近的一个根。 实验实现：

1、算法：

xk?1?xk??f(xk)下山因子从??1开始，逐次将?减半进行试算，直到能使下

f?(xk)降条件f(xk?1)?f(xk)成立为止。再将得到的xk?1循环求得方程根近似值。

2、程序编写代码如下：

function z=f(x) z=2*x^3-5*x-17;

function z=df(x) z=6*x^2-5;

3、运行过程及结果：

实验体会：

牛顿法构造简单，收敛很快，但对初值选取有要求，且牛顿法每次迭代都要计

冶研1201班 s20120273

18811349978

2

北京科技大学 冶金与生态工程学院 冶金工程

算函数的导数，为了克服上述缺点，产生了牛顿下山法。对于较差初值，牛顿下山法优于牛顿法，但牛顿下山法只有线性收敛。

实验二 矩阵的列主元三角分解

一、 实验目的：

学会矩阵的三角分解，并且能够用MATLAB编写相关程序，实现矩阵的三角分解，解方程组。 二、 实验内容：

?1?2??3??4?5??6?7?112345611123451111234111112311111121??x1??7??8??x?1?????2??10?1??x3??????1??x4???13??17?1??x6??????221??x7????28?1????????

用列主元消去法求解方程组（实现PA=LU) 要求输出： （1）计算解X; （2）L,U;

（3）正整型数组IP(i),(i=1,···,n) (记录主行信息)。 三、 实验实现：

1、算法： 列主元三角分解

和普通Dooliitle分解不同，第k步分解时为了避免用绝对值很小的数uii作除数，

引进量Si?aik??limumk i?k,k?1,?,n

m?1k?1若St?maxSi，则将矩阵的第t行与第k行元素互换，再进行正常的

k?i?n冶研1201班 s20120273 18811349978 3

北京科技大学 冶金与生态工程学院 冶金工程

Doolittle分解。

2、程序代码如下：（编写liezhuyuan.m文件）：

function [RA,RB,n,X]=liezhuyuan(A,b)

B=[A b]; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0,

disp('请注意：因为RA~=RB，所以此方程组无解.') return end

if RA==RB if RA==n

disp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1

[Y,j]=max(abs(B(p:n,p))); C=B(p,:); B(p,:)= B(j+p-1,:); B(j+p-1,:)=C; for k=p+1:n

m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);

end end

b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1

X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end

else

disp('请注意：因为RA=RB

（编写lu.m文件）：

function lu(A) [n,m]=size(A); B=zeros(1,n);

L=eye(n);U(1,:)=A(1,:); IP=zeros(n,1); IP(1,1)=U(1,1); for k=2:n

[m,p]=max(abs(A(k:n,k)));%p是m的行标 if m==0;disp('A奇异'),return;end if p>k

B(k:n)=A(k,k:n);A(k,k:n)=A(p,k:n); A(p,k:n)=B(k:n); end

L(k:n,k-1)=A(k:n,k-1)/U(k-1,k-1);

冶研1201班 s20120273

18811349978

4

北京科技大学 冶金与生态工程学院 冶金工程

U(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n); if k

A(k+1:n,k)=A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k); end

IP(k,1)=U(k,k); end disp(IP) disp(L) disp(U)

在Command Windows中输入：

A=[1 1 1 1 1 1 1;2 1 1 1 1 1 1;3 2 1 1 1 1 1;4 3 2 1 1 1 1;5 4 3 2 1 1 1;6 5 4 3 2 1 1;7 6 5 4 3 2 1]; b=[7;8;10;13;17;22;28]; [RA,RB,n,X] =liezhuyuan(A,b) Lu(A)

3、运行结果：

冶研1201班 s20120273

18811349978

5

北京科技大学 冶金与生态工程学院 冶金工程

四．实验体会：通过本次实验，我掌握了消元法解方程的一些基本算法以及用matlab实现矩阵的几种基本计算，对MATLAB软件有了更深的了解。

冶研1201班 s20120273 18811349978 6

北京科技大学 冶金与生态工程学院 冶金工程

实验三 SOR迭代法

实验内容：

?10x1?x2?9?采用SOR方法求解方程??x1?10x2?2x3?7

??2x?10x?623?采用0向量为初始向量，迭代以x(k?1)?x(k)??10?6结束，但迭代次数不操作1000

次，松弛因子使用0.1*i 其中5?i?18，给出使用的各松弛因子及对应迭代次数。 1. 算法：x（k?1)i?（1??）x（k)i??aii(bi??aijxj?1i?1(k?1)j?j?i?1?anijx(jk))

x(k+1)=(D+ωL)-1[（1-ω）D-ωU] x(k)+ ω (D+ωL)-1b

2.程序：编写M函数文件f.m，如下：

function y=f(w)

A=[10 -1 0;-1 10 -2;0 -2 10]; b=[9;7;6]; X0=[0;0;0]; X=X0;

for K=1:1000

for i=1:3

X(i)=w*(b(i)-A(i,:)*X)/A(i,i)+X(i); end

if norm(X-X0)<0.000001 disp(X);disp(K) return; end X0=X; end

3.调用f.m，分别求出迭代结果，如下： >> f(0.5)

0.9958 0.9579 >> f(0.6)

0.9958 0.9579 >> f(0.7)

冶研1201班 s20120273

18811349978

7

0.7916 26

0.7916 21

北京科技大学 冶金与生态工程学院 冶金工程

0.9958 0.9579 0.7916 16

>> f(0.8)

0.9958 0.9579 0.7916 13

>> f(0.9)

0.9958 0.9579 0.7916

10

>> f(1.0)

0.9958 0.9579 0.7916 7

>> f(1.1)

0.9958 0.9579 0.7916 8

>> f(1.2)

0.9958 0.9579 0.7916 11

>> f(1.3)

0.9958 0.9579 0.7916 13

>> f(1.4)

0.9958 0.9579 0.7916 17

>> f(1.5)

0.9958 0.9579 0.7916 23

>> f(1.6)

0.9958 0.9579 0.7916 30

>> f(1.7)

0.9958 0.9579 0.7916 43

>> f(1.8)

0.9958 0.9579 0.7916 68

冶研1201班 s20120273 18811349978

8

北京科技大学 冶金与生态工程学院 冶金工程

实验四 多项式插值

实验内容： 下列数据点的插值

x y 0 0 1 1 4 2 9 3 16 4 25 5 36 6 49 7 64 8 可以得到平方根函数的近似，在区间[0,64]上作图. (注：画图时以步长0.1计算出各点插值，再画出641个点的图，若能给出具体解析式，直接画出连续图)

(1) 用这九个点作8次多项式插值L8(x)，作图，从得到的结果看在[0,64]上，哪个插值区间更精确？

(2) 使用三次插值（不是三弯矩，每个点采用距离最近的四个结点进行二次插值，如插值点为10，则取1,4,9,16四点做三次插值），作图，从得到的结果看在[0,64]上，哪个插值区间更精确？

(3) 比较（1）（2）的结果。

实验实现：

(1).多项式插值程序：

x=[0 1 4 9 16 25 36 49 64]; y=[0 1 2 3 4 5 6 7 8]; xx=0:0.1:64;

yy=interp1(x,y,xx); plot(xx,yy); xlabel('x'); ylabel('y');

冶研1201班 s20120273 18811349978 9

北京科技大学 冶金与生态工程学院 冶金工程

(2).三次插值程序：

x=[0 1 4 9 16 25 36 49 64]; y=[0 1 2 3 4 5 6 7 8]; xx=0:0.1:64; yy=spline(x,y,xx); plot(xx,yy); xlabel('x'); ylabel('y'); grid on 结果：

冶研1201班 s20120273 18811349978 10

北京科技大学 冶金与生态工程学院 冶金工程

(3).比较10开平方后结果为3.1623,则第一种比较精确。 实验体会：

通过本次实验，学会了多项式插值方法，对MATLAB软件有了更深的了解，更加熟练使用matlab程序。

冶研1201班 s20120273 18811349978 11

北京科技大学 冶金与生态工程学院 冶金工程

实验五 数据拟合

实验内容：

1 浓度变化规律

在化学反应中，为研究某化合物的浓度随时间的变化规律，测得一组数据如下：

t(分) 1 2 3 4 5 6 7

y t(分) y

4 9 10

6.4 10 10.2

8.0 11 10.32

8.4 12 10.42

9.28 13 10.5

9.5 14 10.55

9.7 15 10.58

8 9.86 16 10.6

计算出y与t的二次、三次多项式拟合函数，并推断20，40分钟的浓度值，

解：二次多项式拟合：

即最小二乘拟合二次多项式为：y=-0.0445t2+1.0711t+4.3252 再输入：

即f(20)=7.9502,f(40)=-24.0172. 三次多项式拟合：

冶研1201班 s20120273

18811349978

12

北京科技大学 冶金与生态工程学院 冶金工程

即最小二乘拟合三次多项式为：y=0.0060t3-0.1963t2+2.1346t+2.5952 再输入：

即f(20)=14.3973,f(40)=154.8653.

2. 经济增长模型（注：可转变为矛盾方程组求解）

增加生产、发展经济所以靠的主要因素有增加投资、增加劳动力以及技术革新等，在研究国民经济产值与这些因素的数量关系时，由于技术水平不像资金、劳动力那样容易定量化，作为初步的模型，可认为技术水平不变，只讨论产值和资金、劳动力之间的关系。在科学技术发展不快时，如资本主义经济发展的前期，这种模型是有意义的。

用Q，K，L分别表示产值、资金、劳动力，要寻求的数量关系Q（K,L）。经过简化假设与分析，在经济学中，推导出一个著名的Cobb-Douglas生产函数

Q(K,L)?aK?L? （*）

式中?,?,a要经由经济统计α数据确定。现有美国马萨诸塞州1900-1926年上述三个经济指数的统计数据，见下表，试用数据拟合的方法，求出（*）式中的参数?,?,a。

马萨诸塞1900-1926年统计数据 t 1900 1901 1902 1903 Q 1.05 1.18 1.29 1.3 K 1.04 1.06 1.16 1.22 L 1.05 1.08 1.18 1.22 18811349978

13

冶研1201班 s20120273

北京科技大学 冶金与生态工程学院 冶金工程

1904 1905 1906 1907 1908 1909 1910 1911 1912 1913 1914 1915 1916 1917 1918 1919 1920 1921 1922 1923 1924 1925 1926 1.3 1.42 1.5 1.52 1.46 1.6 1.69 1.81 1.93 1.95 2.01 2 2.09 1.96 2.2 2.12 2.16 2.08 2.24 2.56 2.34 2.45 2.58 1.27 1.37 1.44 1.53 1.57 2.05 2.51 2.63 2.74 2.82 3.24 3.24 3.61 4.1 4.36 4.77 4.75 4.54 4.54 4.58 4.58 4.58 4.54 1.17 1.3 1.39 1.47 1.31 1.43 1.58 1.59 1.66 1.68 1.65 1.62 1.86 1.93 1.96 1.95 1.9 1.58 1.67 1.82 1.6 1.61 1.64

解：对函数Q(K,L)?aK?L?取对数得：lnQ?lna??lnK??lnL令

A??1lnKlnL?，x??lna??? lnQ?Ax，从而x?lnQ，即待定参数可用解超A定方程组来求解。

>> Q=[1.05 1.18 1.29 1.3 1.3 1.42 1.5 1.52 1.46 1.6 1.69 1.81 1.93 1.95 2.01 2 2.09 1.96 2.2 2.12 2.16 2.08 2.24 2.56 2.34 2.45 2.58];

K=[1.04 1.06 1.16 1.22 1.27 1.37 1.44 1.53 1.57 2.05 2.51 2.63 2.74 2.82 3.24 3.24 3.61 4.1 4.36 4.77 4.75 4.54 4.54 4.58 4.58 4.58 4.54];

L=[1.05 1.08 1.18 1.22 1.17 1.3 1.39 1.47 1.31 1.43

冶研1201班 s20120273

18811349978

14

北京科技大学 冶金与生态工程学院 冶金工程

1.58 1.59 1.66 1.68 1.65 1.62 1.86 1.93 1.96 1.95 1.9 1.58 1.67 1.82 1.6 1.61 1.64];

>> A=[ones(size(Q)),log(K),log(L)]; >> x=A\\log(Q) x = 0.1626 0.4153 0.0619 >> a=exp(x(1)) a = 1.1766

所以a =1.1766，?=0.4153，?=0.0619

实验体会：

通过本次实验，学会了数据拟合的方法，对MATLAB软件有了更深的了解，更加熟练使用matlab程序。

冶研1201班 s20120273 18811349978 15

北京科技大学 冶金与生态工程学院 冶金工程

实验六 数值积分

实验目的

掌握复化梯形公式和复化Simpson公式，会编写程序，上机进行实例计算。 实验内容

计算：

ln(1?x)?0xdx

1（1）编写复化梯形公式和复化Simpson公式通用子程序，，分别采用4,8,16,32,64等分区间计算。

（2）使用Romberg求积公式。 实验步骤：

1、算法：

（1）复化梯形算法：

（2）复化simpson算法：

（3）龙贝格算法：

2、程序代码如下：

function fuhuatixing(a,b,n) h=(b-a)/n; m=0;

for i=1:n-1

m=log(1+a+i*h)/(a+i*h)+m; end

t=h*(log(1+a)/a+2*m+log(1+b)/b)/2; disp(t)

function simpson(a,b,n) h=(b-a)/(2*n);

m=log(1+a+h)/(a+h);

冶研1201班 s20120273

18811349978

16

北京科技大学 冶金与生态工程学院 冶金工程

r=0;t=0; for i=1:n-1

r=log(1+a+(2*i+1)*h)/(a+(2*i+1)*h)+r; t=log(1+a+2*i*h)/(a+2*i*h)+t; end

s=h*(log(1+a)/a+4*(m+r)+2*t+log(1+b)/b)/3; disp(s);

function [R,k,T]=romberg(fun,a,b,tol) k=0; n=1; h=b-a;

T=h/2*(fun(a)+fun(b)); err=1; while err>=tol k=k+1; h=h/2; tmp=0; for i=1:n

tmp=tmp+fun(a+(2*i-1)*h); end

T(k+1,1)=T(k)/2+h*tmp; for j=1:k

T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1); end n=n*2;

err=abs(T(k+1,k+1)-T(k,k)); end R=T(k+1,4);

3、运行结果：

冶研1201班 s20120273 18811349978 17

北京科技大学 冶金与生态工程学院 冶金工程

冶研1201班 s20120273 18811349978 18

北京科技大学 冶金与生态工程学院 冶金工程

实验体会：

通过本次实验，学会了用复化梯形求积公式和simpson公式，以及龙贝格求积公式求解数值积分的方法，对MATLAB软件有了更深的了解，更加熟练的使用matlab程序。

冶研1201班 s20120273 18811349978 19

北京科技大学 冶金与生态工程学院 冶金工程

实验七 数值积分

实验内容

利用改进Euler法，经典四级四阶Rung-Kutta方法求其数值解，：

??4x?xy?y? ?y??y(0)?1 分别取h=0.2,0.1,0.05时数值解。 1. 改进Euler法 程序：

for n=1:10 x=0;y=1;h=0.2; for i=1:n

y1=y+h*(4*x/y-x*y);

y2=y+h*(4*(x+h)/(y1)-(x+h)*(y1)); y=(y1+y2)/2; x=x+h; end p(n)=y; end p

h=0.2

冶研1201班 s20120273 18811349978 0?x?2

y?4?3e?x2）20

（注：初值问题的精确解

北京科技大学 冶金与生态工程学院 冶金工程

h=0.1

h=0.05

冶研1201班 s20120273 18811349978

21

北京科技大学 冶金与生态工程学院 冶金工程

2. RK法： 程序：

for n=1:10 x=0;y=1;h=0.2; for i=1:n

k1=h*(4*x/y-x*y);

k2=h*(4*(x+h/2)/(y+k1/2)-(x+h/2)*(y+k1/2)); k3=h*(4*(x+h/2)/(y+k2/2)-(x+h/2)*(y+k2/2)); k4=h*(4*(x+h)/(y+k3/2)-(x+h)*(y+k3)); y=y+(k1+2*k2+2*k3+k4)/6; x=x+h; end p(n)=y; end p

h=0.2

冶研1201班 s20120273 18811349978 22

北京科技大学 冶金与生态工程学院 冶金工程

h=0.1

h=0.05

冶研1201班 s20120273 18811349978

23

北京科技大学 冶金与生态工程学院 冶金工程

实验体会：

通过本次实验，学会了用改进欧拉方法和经典R-K方法求解常微分方程数值解的方法，对MATLAB软件有了更深的了解，更加熟练使用matlab程序。

冶研1201班 s20120273 18811349978 24

本文来源：https://www.bwwdw.com/article/wr3x.html