计算方法大作业共21页.docx
《计算方法大作业共21页.docx》由会员分享，可在线阅读，更多相关《计算方法大作业共21页.docx（25页珍藏版）》请在冰豆网上搜索。
计算方法大作业共21页
实验（shíyàn）一牛顿（niúdùn）下山法
一、实验（shíyàn）目的：
1、掌握牛顿下山法求解（qiújiě）方程根的推导原理。
2、理解牛顿下山法的具体算法与相应程序的编写。
二、实验内容：
采用牛顿下山法求方程2x3-5x-17=0在2附近的一个根。
三、实验实现：
1、算法：
下山因子从
开始，逐次将
减半进行试算，直到能使下降条件
成立为止。
再将得到的
循环求得方程根近似值。
2、程序代码如下：
function[p,k]=NewtonDownHill（f,df,p0）
N=2000;Tol=10^（-5）;e=10^（-8）;
fork=1:
N
lamda=1;
p1=p0-lamda*f（p0）/df（p0）;
while（abs（f（p1））>=abs（f（p0））&lamda>e）
lamda=lamda/2;
p1=p0-lamda*f（p0）/df（p0）;
end
ifabs（p1-p0）end
p0=p1;
end
ans=p1
3、运行（yùnxíng）结果：
四、实验（shíyàn）体会：
牛顿下山法可以较快求的方程（fāngchéng）结果，对于该题，只需要5步。
运用计算机的数值迭代法可以很快求得满足精度要求的结果。
实验（shíyàn）二矩阵的列主元三角分解
一、实验目的：
学会矩阵的三角分解，并且能够用MATLAB编写相关程序，实现矩阵的三角分解，解方程组。
二、实验内容：
用列主元消去法求解方程组（实现PA=LU）要求输出：
（1）计算解X;
（2）L,U;
（3）正整型数组IP（i）,（i=1,···,n）（记录（jìlù）主行信息）。
三、实验（shíyàn）实现：
1、算法（suànfǎ）：
列主元三角分解和普通Dooliitle分解不同（bùtónɡ），第k步分解时为了避免用绝对值很小的数
作除数，
引进量
若
，则将矩阵的第t行与第k行元素互换，再进行正常的Doolittle分解。
2、程序代码如下：
clearall;clc;
A=[1111111;
2111111;
3211111;
4321111;
5432111;
6543211;
7654321];
b=[7;8;10;13;17;22;28];
n=length（A）;
IP=eye（n）;
U=zeros（n）;
L=eye（n）;
[m,p]=max（A（:
1））;
C1=b
（1）;b
（1）=b（p）;b（p）=C1;
C2（1:
n）=IP（1,:
）;IP（1,:
）=IP（p,:
）;IP（p,:
）=C2（1:
n）;
C2（1:
n）=A（1,:
）;A（1,:
）=A（p,:
）;A（p,:
）=C2（1:
n）;
U（1,:
）=A（1,:
）;
fork=2:
n
L（k:
n,k-1）=A（k:
n,k-1）/U（k-1,k-1）;
j=1;
fori=k:
n
S（j）=A（i,k）-L（i,1:
k-1）*U（1:
k-1,k）;
j=j+1;
end
[m,p]=max（abs（S（:
）））;
C2（1:
n）=A（k,1:
n）;A（k,1:
n）=A（k+p-1,1:
n）;A（k+p-1,1:
7）=C2（1:
n）;
C2（1:
n）=IP（k,1:
n）;IP（k,1:
n）=IP（k+p-1,1:
n）;IP（k+p-1,1:
n）=C2（1:
n）;
U（k,k:
n）=A（k,k:
n）-L（k,1:
k-1）*U（1:
k-1,k:
n）;
ifkA（k+1:
n,k）=A（k+1:
n,k）-L（k+1:
n,1:
k-1）*U（1:
k-1,k）;
end
end
disp（'L='）;disp（L）;
disp（'U='）;disp（U）;
disp（'IP='）;disp（IP）;
y=L\b;
x=U\y;
disp（'y='）;disp（y）;
disp（'x='）;disp（x）;
4、运行（yùnxíng）结果：
四、实验（shíyàn）体会：
列主元三角分解（fēnjiě）原理很简单，用MATLAB实现起来却不容易。
对MATLAB还需要深入学习。
实验（shíyàn）三SOR迭代法
一、实验（shíyàn）内容：
采用SOR方法求解方程
采用0向量为初始向量，迭代以
结束，但迭代次数不操作1000次，松弛因子使用
其中
，给出使用的各松弛因子及对应迭代次数。
二、实验实现：
1、程序代码如下：
A=[10-10;
-110-2;
0-210];
b=[9;7;6];
x0=[0;0;0];e=10^（-6）;n=1000;
x=x0;
j=18;
w=0.1*j;
fork=1:
n
fori=1:
length（b）
x（i）=w*（b（i）-A（i,:
）*x）/A（i,i）+x（i）;
ifnorm（x-x0）disp（w）;disp（x'）;disp（k）;disp（'------------'）
return
end
x0=x;
end
end
2、运行结果（jiēguǒ）如下：
上图显示（xiǎnshì）为松弛因子分别为0.5、1、1.5和1.8时的结果（jiēguǒ）。
实验（shíyàn）四多项式插值
一、实验（shíyàn）内容：
下列（xiàliè）数据点的插值
x
0
1
4
9
16
25
36
49
64
y
0
1
2
3
4
5
6
7
8
可以得到平方根函数的近似，在区间[0,64]上作图.（注：
画图时以步长0.1计算（jìsuàn）出各点插值，再画出641个点的图，若能给出具体解析式，直接画出连续图）
（1）用这九个点作8次多项式插值L8（x）.，作图，从得到的结果看在[0,64]上，哪个插值区间更精确？
（2）使用三次插值（不是三弯矩，每个点采用距离最近的四个结点进行三次插值，如插值点为10，则取1,4,9,16四点做三次插值），作图，从得到的结果看在[0,64]上，哪个插值区间更精确？
（3）比较
（1）
（2）的结果。
二、实验实现：
1、程序代码如下：
x=[01491625364964];
y=[012345678];
x1=0:
0.1:
64;
y1=zeros（641,1）;y2=zeros（641,1）;
forj=1:
1:
641
fork=1:
9
t=1;
fori=1:
9
ifi~=k
t=t*（x1（j）-x（i））/（x（k）-x（i））;
end
end
y1（j）=y1（j）+t*y（k）;
end
end
x2=0:
0.1:
64;
forj=1:
41
fork=1:
4
t=1
fori=1:
4
ifi~=k
t=t*（x2（j）-x（i））/（x（k）-x（i））;
end
end
y2（j）=y2（j）+t*y（k）;
end
end
forj=42:
1:
91
fork=2:
5
t=1
fori=2:
5
ifi~=k
t=t*（x2（j）-x（i））/（x（k）-x（i））;
end
end
y2（j）=y2（j）+t*y（k）;
end
end
forj=92:
1:
161
fork=3:
6
t=1
fori=3:
6
ifi~=k
t=t*（x2（j）-x（i））/（x（k）-x（i））;
end
end
y2（j）=y2（j）+t*y（k）;
end
end
forj=162:
1:
251
fork=4:
7
t=1
fori=4:
7
ifi~=k
t=t*（x2（j）-x（i））/（x（k）-x（i））;
end
end
y2（j）=y2（j）+t*y（k）;
end
end
forj=252:
1:
361
fork=5:
8
t=1
fori=5:
8
ifi~=k
t=t*（x2（j）-x（i））/（x（k）-x（i））;
end
end
y2（j）=y2（j）+t*y（k）;
end
end
forj=362:
1:
641
fork=6:
9
t=1
fori=6:
9
ifi~=k
t=t*（x2（j）-x（i））/（x（k）-x（i））;
end
end
y2（j）=y2（j）+t*y（k）;
end
end
subplot（2,1,1）
plot（x,y,'ro',x1,y1,'b.',x1,sqrt（x1）,'r'）
gridon;axis（[070010]）;
subplot（2,1,2）
plot（x,y,'ro',x2,y2,'k.',x1,sqrt（x1）,'r'）
gridon;axis（[070010]）;
2、运行（yùnxíng）结果如下（rúxià）：
3、结果（jiēguǒ）分析
（1）由九个点做八次（bācì）插值，由上图可知，仅在区间[016]上较精确，其他区间均与真实值有较大差距，区间[2564]上，甚至出现严重波动。
（2）做分段三次插值，发现在各个区间上，都能获得比较精确（jīngquè）的插值。
仅在区间[04]上有少许偏差。
（3）两个结果对比说明，由于龙格现象的存在，高次插值虽能在更多点上与被插值函数的函数值相同，但从整体（zhěngtǐ）上看不一定能取得很好的逼近效果，甚至出现震荡。
这时选用分段低次插值可以取得较好的效果。
实验五数据拟合
一、实验内容：
1、浓度变化规律
在化学反应中，为研究某化合物的浓度随时间的变化规律，测得一组数据如下：
t（分）
1
2
3
4
5
6
7
8
y
4
6.4
8.0
8.4
9.28
9.5
9.7
9.86
t（分）
9
10
11
12
13
14
15
16
y
10
10.2
10.32
10.42
10.5
10.55
10.58
10.6
计算出y与t的二次、三次多项式拟合（nǐhé）函数，并推断20，40分钟的浓度值，
（1）程序代码如下（rúxià）：
t=1:
1:
16;
y=[46.48.08.49.289.59.79.861010.210.3210.4210.510.5510.5810.6];
x=1:
0.1:
50;
y2=polyfit（t,y,2）;
y3=polyfit（t,y,3）;
y_2=poly2str（y2,'t'）
y_3=poly2str（y3,'t'）
y22=polyval（y2,x）;
y33=polyval（y3,x）;
y2_20=polyval（y2,20）
y2_40=polyval（y2,40）
y3_20=polyval（y3,20）
y3_40=polyval（y3,40）
plot（x,y22,'r',x,y33,'b',t,y,'ko'）;holdon
plot（20,y2_20,'b*',40,y2_40,'b*'）;holdon
plot（20,y3_20,'b*',40,y3_40,'b*'）;holdon
gtext（'二次拟合（nǐhé）'）;
gtext（'三次（sāncì）拟合'）;
（2）运行结果如下：
二次、三次拟合多项式为：
二次、三次拟合情况下，20分钟和40分钟的浓度值为：
二次、三次（sāncì）拟合的拟合曲线为：
2、经济增长模型（注：
可转变为矛盾（máodùn）方程组求解）
增加生产、发展经济所以靠的主要因素有增加投资、增加劳动力以及技术革新等，在研究国民经济产值与这些因素的数量关系时，由于技术水平不像资金、劳动力那样容易（róngyì）定量化，作为初步的模型，可认为技术水平不变，只讨论产值和资金、劳动力之间的关系。
在科学技术发展不快时，如资本主义经济发展的前期，这种模型是有意义的。
用Q，K，L分别表示产值（chǎnzhí）、资金、劳动力，要寻求的数量关系Q（K,L）。
经过简化假设与分析，在经济学中，推导出一个著名的Cobb-Douglas生产函数
（*）
式中
要经由经济（jīngjì）统计α数据确定。
现有美国马萨诸塞州1900-1926年上述三个经济指数的统计数据，见下表，试用数据拟合的方法（fāngfǎ），求出（*）式中的参数
。
马萨诸塞1900-1926年统计数据
t
Q
K
L
1900
1.05
1.04
1.05
1901
1.18
1.06
1.08
1902
1.29
1.16
1.18
1903
1.3
1.22
1.22
1904
1.3
1.27
1.17
1905
1.42
1.37
1.3
1906
1.5
1.44
1.39
1907
1.52
1.53
1.47
1908
1.46
1.57
1.31
1909
1.6
2.05
1.43
1910
1.69
2.51
1.58
1911
1.81
2.63
1.59
1912
1.93
2.74
1.66
1913
1.95
2.82
1.68
1914
2.01
3.24
1.65
1915
2
3.24
1.62
1916
2.09
3.61
1.86
1917
1.96
4.1
1.93
1918
2.2
4.36
1.96
1919
2.12
4.77
1.95
1920
2.16
4.75
1.9
1921
2.08
4.54
1.58
1922
2.24
4.54
1.67
1923
2.56
4.58
1.82
1924
2.34
4.58
1.6
1925
2.45
4.58
1.61
1926
2.58
4.54
1.64
（1）程序代码
对（*）式两边取对数（duìshù）并整理得：
令A=[1lnKlnL],
则lnQ=Ax,x=A\lnQ
K=[1.041.061.161.221.271.371.441.531.572.052.512.632.742.823.243.243.614.14.364.774.754.544.544.584.584.584.54]';
L=[1.051.081.181.221.171.31.391.471.311.431.581.591.661.681.651.621.861.931.961.951.91.581.671.821.61.611.64]';
Q=[1.051.181.291.31.31.421.51.521.461.61.691.811.931.952.0122.091.962.22.122.162.082.242.562.342.452.58]';
lnK=log（K）;lnL=log（L）;lnQ=log（Q）;
I=ones（size（Q））;
A=[IlnKlnL];
x=A\lnQ
a=exp（x
（1））
alpha=x
（2）
beta=x（3）
（2）运行（yùnxíng）结果如下：
实验（shíyàn）六数值积分
一、实验（shíyàn）目的
掌握（zhǎngwò）复化梯形公式和复化Simpson公式，会编写程序，上机进行（jìnxíng）实例计算。
二、实验内容
计算
（1）编写复化梯形公式和复化Simpson公式通用子程序，，分别采用4,8,16,32,64等分区间计算。
程序代码如下：
formatlong
f=inline（'log（1+x）./x'）;
a=10^-8;b=1;
fori=2:
6
n=2^i;
h=（b-a）/n;
x=zeros（n-1）;sum=0;
fori=1:
1:
n-1
x（i）=a+i*h;
sum=sum+f（x（i））;
end
fa=f（a）;fb=f（b）;
T=h/2*（fa+2*sum+fb）;
n;
T
end
运行（yùnxíng）结果为：
（2）使用（shǐyòng）Romberg求积公式。
程序代码如下（rúxià）：
f=inline（'log（1+x）./x'）;
a=10^-8;b=1;
M=1;h=b-a;J=0;err=1;R=zeros（4,4）;eps=10^-8;
R（1,1）=h*（f（a）+f（b））/2;
while（err>eps）
J=J+1;
h=h/2;
x=a+h:
2*h:
b-h;
R（J+1,1）=R（J,1）/2+h*sum（f（x））;
forK=1:
min（3,J）
R（J+1,K+1）=R（J+1,K）+（R（J+1,K）-R（J,K））/（4^K-1）;
end
if（J>3）
err=abs（R（J+1,4）-R（J,4））;
end
end
S=R（J,4）
运行（yùnxíng）结果为：
实验（shíyàn）七数值积分
实验（shíyàn）内容
利用（lìyòng）改进Euler法，经典四级四阶Rung-Kutta方法求其数值解，：
分别取h=0.2,0.1,0.05时数值解。
（注：
初值问题的精确解
）
1、改进Euler法
（1）程序代码如下：
x0=0;y0=1;xn=2;
f=inline（'4*x./y-x.*y'）;
H=[0.20.10.05];
forj=1:
3
h=H（j）;
n=（xn-x0）/h;
x=x0;y=y0;
fori=1:
n
y1=y+h*f（x,y）;
y2=y+h*f（x+h,y1）;
y=（y1+y2）/2;
x=x+h;
Y（j,i）=y;
end
h
Y（j,1:
n）
disp（'_____________'）;
（2）运行结果（jiēguǒ）如下：
2、经典（jīngdiǎn）R-K方法
（1）程序代码如下（rúxià）：
x0=0;y0=1;xn=2;
f=inline（'4*x./y-x.*y'）;
H=[0.20.10.05];
forj=1:
3
h=H（j）;
n=（xn-x0）/h;
x=x0;y=y0;
fori=1:
n
K1=h*f（x,y）;
K2=h*f（x+h/2,y+K1/2）;
K3=h*f（x+h/2,y+K2/2）;
K4=h*f（x+h,y+K3）;
y=y+（K1+2*K2+2*K3+K4）/6;
x=x+h;
Y（j,i）=y;
end
h
Y（j,1:
n）
disp（'_____________'）;
end
（2）运行结果（jiēguǒ）如下：
内容摘要
（1）实验一牛顿下山法
实验目的：
掌握牛顿下山法求解方程根的推导原理
（2）IP（p,:
）=C2（1:
n）
（3）I=ones（size（Q））
（4）x=A\lnQ
a=exp（x
（1））
alpha=x
（2）
beta=x（3）
（2）运行结果如下：
实验六数值积分
实验目的
掌握复化梯形公式和复化Simpson公式，会编写程序，上机进行实例计算
（5）程序代码如下：
formatlong
f=inline（'log（1+x）./x'）
（6）sum=sum+f（x（i））
（7）程序代码如下：
f=inline（'log（1+x）./x'）
（8）x=x+h
（9）Y（j,i）=y
升级会员 