文档视界 最新最全的文档下载
当前位置:文档视界 › 数值计算方法作业

数值计算方法作业

数值计算方法作业
数值计算方法作业

目录

1非线性方程求根 (1)

1.1迭代法 (1)

1.2牛顿法 (1)

1.3弦截法 (2)

1.4二分法 (3)

2插值 (4)

2.1线性插值 (4)

2.2二次插值 (5)

2.3拉格朗日插值 (6)

2.4分段线性插值 (7)

2.5分段二次插值 (8)

3数值积分 (10)

3.1复化矩形积分法 (10)

3.2 复化梯形积分法 (10)

3.3 复化辛甫生积分法 (11)

3.4 变步长梯形积分法 (12)

4线性方程组数值解法 (13)

4.1约当消去法 (13)

4.2高斯消去法 (14)

4.3 三角分解法 (15)

4.4 雅克比迭代法 (17)

4.5高斯-塞德尔迭代法 (19)

5常微分方程数值解 (20)

5.1 显式欧拉公式 (20)

5.2 欧拉公式预测校正系统 (21)

5.3 两步欧拉公式 (22)

5.4 改进欧拉公式 (23)

5.5 四阶龙格-库塔法 (24)

1非线性方程求根

1.1迭代法

1.1.1程序代码

Private Sub Command1_Click()

x0 = Val(InputBox("请输入初始值X0:"))

ep = Val(InputBox("请输入误差限EP:"))

s = 0

While s = 0

x = 2 * Sin(x0)

If Abs(x - x0) < ep Then

Print "x=", x

s = 1

Else

x0 = x

End If

Wend

End Sub

1.1.2 算例

求方程2sin(x)-x=0的根,误差限为0.0001

1.1.3运行结果

1.2牛顿法

1.2.1程序代码

Private Sub Command1_Click()

x0 = Val(InputBox("请输入初始值X0:"))

ep = Val(InputBox("请输入误差限EP:"))

S = 0

While S = 0

X1 = x0 - (x0 ^ 2 * Exp(x0) - 8 * x0) / (x0 ^ 2 * Exp(x0) + 2 * x0 * Exp(x0) - 8) If Abs(X1 - x0) < ep Then

Print "非线性方程的根X="; X1

S = 1

Else

x0 = X1

End If

Wend

End Sub

1.2.2算例

求方程的最小正根,误差限为0.0001 1.2.3运行结果

1.3弦截法

1.3.1程序代码

Private Sub Command1_Click()

Dim x0, x1, n, ep, y0, y1

x0 = Val(InputBox("请输入区间左端点:"))

x1 = Val(InputBox("请输入区间右端点:"))

n = Val(InputBox("最大迭代数n:"))

ep = Val(InputBox("请输入误差限:"))

s = 0

f = 0

While f = 0

y0 = x0 ^ 3 - 3 * x0 ^ 2 + 2 * x0 - 18

y1 = x1 ^ 3 - 3 * x1 ^ 2 + 2 * x0 - 18

x = x1 - y1 * (x1 - x0) / (y1 - y0)

If Abs(x - x1) < ep Then

Print "x="; x

f = 1

Else

s = s + 1

If s > n Then

Print "计算失败"

f = 1

Else

x0 = x1

x1 = x

End If

End If

Wend

End Sub

1.3.2 算例

求x3-3x2+2x-18=0在[3,5]的根,最大迭代数1000,误差限为0.0001 1.3.3运行结果

1.4二分法

1.4.1程序代码

Private Sub Command1_Click()

a = Val(InputBox("请输入区间左端点a:"))

b = Val(InputBox("请输入区间右端点b:"))

ep = Val(InputBox("请输入误差限EP:"))

s = 0

While s = 0

x = (a + b) / 2

fx = x ^ 3 - 2 * x ^ 2 + 3 * x - 12

fa = a ^ 3 - 2 * a ^ 2 + 3 * a - 12

If fx = 0 Then

Print "x=", x

s = 1

Else

If fx * fa > 0 Then

a = x

Else

b = x

End If

If Abs(b - a) < ep Then

x = (a + b) / 2

Print "x="; x

s = 1

Else

End If

End If

Wend

End Sub

1.4.2 算例

求x3-2x2+3x-12=0在[1,3]的根,误差限为0.0001 1.4.3 运行结果

2插值

2.1线性插值

2.1.1 程序代码

Private Sub Command1_Click()

X0 = Val(InputBox("请输入插值结点X0:"))

y0 = Val(InputBox("请输入插值结点Y0:"))

X1 = Val(InputBox("请输入插值结点X1:"))

Y1 = Val(InputBox("请输入插值结点Y1:"))

f = 0

While f = 0

x = Val(InputBox("请输入插值点X:"))

L0 = (x - X1) / (X0 - X1)

L1 = (x - X0) / (X1 - X0)

y = L0 * y0 + L1 * Y1

Print "x="; x, "y="; y

f = InputBox("是否继续(0/1)?")

Wend

End Sub

2.1.2算例

已知(16,5)、(36,7),求18处的值

2.1.3 运行结果

2.2二次插值

2.2.1 程序代码

Private Sub Command1_Click()

X0 = Val(InputBox("请输入插值结点X0:"))

y0 = Val(InputBox("请输入插值结点Y0:"))

X1 = Val(InputBox("请输入插值结点X1:"))

Y1 = Val(InputBox("请输入插值结点Y1:"))

X2 = Val(InputBox("请输入插值结点X2:"))

Y2 = Val(InputBox("请输入插值结点Y2:"))

f = 0

While f = 0

x = Val(InputBox("请输入插值点X:"))

L0 = (x - X1) * (x - X2) / ((X0 - X1) * (X0 - X2))

L1 = (x - X0) * (x - X2) / ((X1 - X0) * (X1 - X2))

L2 = (x - X0) * (x - X1) / ((X2 - X0) * (X2 - X1))

y = L0 * y0 + L1 * Y1 + L2 * Y2

Print "x="; x, "y="; y

f = InputBox("是否继续(0/1)?")

Wend

End Sub

2.2.2 算例

已知(27,4)、(125,6)、(343,8),求50处的值2.2.3 运行结果

2.3拉格朗日插值

2.3.1 程序代码

Private Sub Command1_Click()

Dim x(), y()

n = Val(InputBox("请输入插值结点数N:"))

ReDim x(n), y(n)

For i = 0 To n

x(i) = Val(InputBox("请输入插值结点X(" + Str(i) + "):"))

y(i) = Val(InputBox("请输入插值结点Y(" + Str(i) + "):"))

Next i

f = 0

While f = 0

xx = Val(InputBox("请输入插值点X:"))

Sum = 0

For i = 0 To n

T = 1

For j = 0 To n

If i <> j Then

T = T * (xx - x(j)) / (x(i) - x(j))

Else

End If

Next j

Sum = Sum + T * y(i)

Next i

Print "x="; xx, "y="; Sum

f = InputBox("是否继续(0/1)?")

Wend

End Sub

2.3.2 算例

已知(27,4)(64,5)(125,6)(216,7)(512,9),N=4,求150处的值

2.3.3 运行结果

2.4分段线性插值

2.4.1 程序代码

Private Sub Command1_Click()

Dim x(), y()

n = Val(InputBox("请输入插值结点数N:"))

ReDim x(n), y(n)

For i = 0 To n

x(i) = Val(InputBox("请输入插值结点X(" + Str(i) + "):"))

y(i) = Val(InputBox("请输入插值结点Y(" + Str(i) + "):"))

Next i

f = 0

While f = 0

xx = Val(InputBox("请输入插值点X:"))

i = 1

a = 0

While a = 0

If xx < x(i) Then

k = i - 1

a = 1

Else

i = i + 1

If i > n - 1 Then

k = n - 1

a = 1

Else

End If

End If

Wend

L0 = (xx - x(k + 1)) / (x(k) - x(k + 1))

L1 = (xx - x(k)) / (x(k + 1) - x(k))

YY = L0 * y(k) + L1 * y(k + 1)

Print "x="; xx, "y="; YY

f = InputBox("是否继续(0/1)?")

Wend

End Sub

2.4.2 算例

已知(256,16)(289,17)(324,18)(361,19),N=3,求315处的值2.4.3 运行结果

2.5分段二次插值

2.5.1 程序代码

Private Sub Command1_Click()

Dim x(), y()

n = Val(InputBox("请输入插值结点数N;"))

ReDim x(n), y(n)

For i = 0 To n

x(i) = Val(InputBox("请输入插值结点X(" + Str(i) + "):"))

y(i) = Val(InputBox("请输入插值结点Y(" + Str(i) + "):"))

Next i

f = 0

While f = 0

xx = Val(InputBox("请输入插值点X:"))

If xx < x(1) Then

k = 1

Else

i = 2

s = 0

While s = 0

If xx < x(i) Then

If xx - x(i - 1) < x(i) - xx Then

k = i - 1

s = 1

Else

k = i

s = 1

End If

Else

i = i + 1

If i > n - 1 Then

k = n - 1

s = 1

Else

End If

End If

Wend

End If

L0 = (xx - x(k)) * (xx - x(k + 1)) / ((x(k - 1) - x(k)) * (x(k - 1) - x(k + 1)))

L1 = (xx - x(k - 1)) * (xx - x(k + 1)) / ((x(k) - x(k - 1)) * (x(k) - x(k + 1)))

L2 = (xx - x(k - 1)) * (xx - x(k)) / ((x(k + 1) - x(k - 1)) * (x(k + 1) - x(k)))

yy = L0 * y(k - 1) + L1 * y(k) + L2 * y(k + 1)

Print "x="; xx, "y="; yy

f = InputBox("是否继续(0/1)?")

Wend

End Sub

2.5.2 算例

已知(196,14)、(169,13)、(144,12)、(121,11),N=3,求150处的值2.5.3 运行结果

3数值积分

3.1复化矩形积分法

3.1.1 程序代码

Private Sub Command1_Click()

a = Val(InputBox("请输入积分下限A:"))

b = Val(InputBox("请输入积分上限B:"))

n = Val(InputBox("请输入积分区间等分数N:"))

h = (b - a) / n

Sum = 0

For i = 1 To n

Sum =Sum + (2 *(a + (i - 1 / 2) * h)) / ((a + (i - 1 / 2) * h) ^ 2 + 1)

Next i

R = h * Sum

Print "复化矩形积分法计算结果R"; R

End Sub

3.1.2 算例

计算等分1000

3.1.3 运行结果

3.2 复化梯形积分法

3.2.1 程序代码

Private Sub Command1_Click()

a = Val(InputBox("请输入积分下限A:"))

b = Val(InputBox("请输入积分上限B:"))

n = Val(InputBox("请输入积分区间等分数N:"))

h = (b - a) / n

Sum = 0

For i = 1 To n - 1

Sum = Sum + (a + i * h) ^ 2 / ((a + i * h) ^ 2 + 1))

Next i

T = h * (a ^ 2 / (a ^ 2 + 1) + b ^ 2 / (b ^ 2 + 1)) / 2 + h * Sum

Print "复化梯形积分法计算结果T:"; T

End Sub

3.2.2 算例

计算等分1000

3.2.3 运行结果

3.3 复化辛甫生积分法

3.3.1 程序代码

Private Sub Command1_Click()

a = Val(InputBox("请输入积分下限A:"))

b = Val(InputBox("请输入积分上限B:"))

n = Val(InputBox("请输入积分区间等分数N:"))

h = (b - a) / n

Sum = 0

For i = 1 To n

Sum = Sum + (a + (i - 1) * h) ^ 4 + 4 * (a + (i - 1 / 2) * h) ^ 4 + (a + i * h) ^ 4 Next i

S = h * Sum / 6

Print "辛普生积分法计算结果S:"; S

End Sub

3.3.2 算例

计算dx 等分100

3.3.3 运行结果

3.4 变步长梯形积分法

3.4.1 程序代码

Private Sub Command1_Click()

a= Val(InputBox("请输入积分下限A:"))

b = Val(InputBox("请输入积分上限B:"))

ep = Val(InputBox("请输入误差限EP:"))

n = 1

h = b - a

T1 = h * (2 * a / (a ^ 3 + 1) + 2 * b / (b ^ 3 + 1)) / 2

f = 0

While f = 0

Sum = 0

For i = 1 To n

Sum = Sum +2 * (a + (i - 1 / 2) * h) / ((a + (i - 1 / 2) * h) ^ 3 + 1)

Next i

T2 = T1 / 2 + h * Sum / 2

If Abs(T2 - T1) < ep Then

Print "变步长梯形积分法计算结果T:"; T2

f = 1

Else

T1 = T2

h = h / 2

n = 2 * n

End If

Wend

End Sub

3.4.2 算例

计算误差限ep=0.0001

3.4.3 运行结果

4线性方程组数值解法

4.1约当消去法

4.1.1 程序代码

Private Sub Command1_Click()

Dim a()

n = Val(InputBox("请输入方程个数N:"))

ReDim a(n, n + 1)

For i = 1 To n

For j = 1 To n + 1

a(i, j) = Val(InputBox("请输入增广矩阵A(" + Str(i) + "," + Str(j) + "):"))

Next j

Next i

For k = 1 To n

m = a(k, k)

For j = k To n + 1

a(k, j) = a(k, j) / m

Next j

For i = 1 To n

If i <> k Then

m = a(i, k)

For j = k To n + 1

a(i, j) = a(i, j) - a(k, j) * m

Next j

End If

Next i

Next k

For i = 1 To n

Print "x(" + Str(i) + ")="; a(i, n + 1)

Next i

End Sub

4.1.2 算例

计算方程组

4.1.3 运行结果

4.2高斯消去法

4.2.1 程序代码

Private Sub Command1_Click()

Dim a(), x()

n = Val(InputBox("请输入方程个数N:"))

ReDim a(n, n + 1), x(n)

For i = 1 To n

For j = 1 To n + 1

a(i, j) = Val(InputBox("请输入増广矩阵A(" + Str(i) + "," + Str(j) + "):")) Next j

Next i

For k = 1 To n - 1

m = a(k, k)

For j = k To n + 1

a(k, j) = a(k, j) / m

Next j

For i = k + 1 To n

m = a(i, k)

For j = k To n + 1

a(i, j) = a(i, j) - a(k, j) * m

Next j

Next i

Next k

x(n) = a(n, n + 1) / a(n, n)

For i = n - 1 To 1 Step -1

Sum = 0

For j = i + 1 To n

Sum = Sum + a(i, j) * x(j)

Next j

x(i) = a(i, n + 1) - Sum

Next i

For i = 1 To n

Print "x(" + Str(i) + ")="; x(i)

Next i

End Sub

4.2.2 算例

计算方程组

4.2.3 运行结果

4.3 三角分解法

4.3.1 程序代码

Private Sub Command1_Click()

Dim a(), b(), L(), U(), x(), Y()

n = Val(InputBox("请输入方程个数N:"))

ReDim a(n, n), L(n, n), U(n, n), b(n), x(n), Y(n)

For i = 1 To n

For j = 1 To n

a(i, j) = Val(InputBox("请输入系数矩阵A(" + Str(i) + "," + Str(j) + "):"))

Next j

b(i) = Val(InputBox("请输入右端常数项B(" + Str(i) + "):"))

Next i

For i = 1 To n

For j = 1 To n

If i > j Then

Sum = 0

For k = 1 To j - 1

Sum = Sum + L(i, k) * U(k, j) Next k

L(i, j) = (a(i, j) - Sum) / U(j, j) Else

Sum = 0

For k = 1 To i - 1

Sum = Sum + L(i, k) * U(k, j) Next k

U(i, j) = a(i, j) - Sum

End If

Next j

Next i

For i = 1 To n

Sum = 0

For j = 1 To i - 1

Sum = Sum + L(i, j) * Y(j) Next j

Y(i) = b(i) - Sum

Next i

For i = n To 1 Step -1

Sum = 0

For j = i + 1 To n

Sum = Sum + U(i, j) * x(j) Next j

x(i) = (Y(i) - Sum) / U(i, i) Next i

For i = 1 To n

Print "x(" + Str(i) + ")="; x(i)

Next i

End Sub

4.3.2 算例

计算方程组

4.3.3 运行结果

4.4 雅克比迭代法

4.4.1 程序代码

Private Sub Command1_Click()

Dim a(), b(), x(), x0()

n = Val(InputBox("方程个数N:"))

ReDim a(n, n), x(n), x0(n), b(n)

Nmax = Val(InputBox("请输入最大迭代次数Nmax:"))

ep = Val(InputBox("请输入误差限EP:"))

For i = 1 To n

For j = 1 To n

a(i, j) = Val(InputBox("请输入系数矩阵A(" + Str(i) + "," + Str(j) + "):"))

Next j

b(i) = Val(InputBox("请输入右端常数项B(" + Str(i) + "):"))

Next i

For i = 1 To n

x0(i) = Val(InputBox("请输入初始值X0(" + Str(i) + "):"))

Next i

For i = 1 To n

For j = 1 To n

Print a(i, j);

Next j

Print b(i)

Next i

m = 0

f = 0

While f = 0

Max = 0

For i = 1 To n

Sum = 0

For j = 1 To n

Sum = Sum + a(i, j) * x0(j) Next j

d = (b(i) - Sum) / a(i, i)

If Abs(d) > Max Then

Max = Abs(d)

Else

End If

x(i) = x0(i) + d

Next i

m = m + 1

If Max < ep Then

For i = 1 To n

Print "x(" + Str(i) + ")="; x(i) Next i

f = 1

Else

If m > Nmax Then

Print "迭代失败"

f = 1

Else

For i = 1 To n

x0(i) = x(i)

Next i

End If

End If

Wend

End Sub

4.4.2 算例

计算方程组ep=0.0001,最大迭数1000,初始值均为0 4.4.3 运行结果

4.5高斯-塞德尔迭代法

4.5.1 程序代码

Private Sub Command1_Click()

Dim a(), b(), x()

n = Val(InputBox("请输入方程个数N:"))

ReDim a(n, n), b(n), x(n)

ep = 0.000001

Nmax = 100000

For i = 1 To n

For j = 1 To n

a(i, j) = Val(InputBox("请输入系数矩阵A(" + Str(i) + "," + Str(j) + "):"))

Next j

b(i) = Val(InputBox("请输入右端常数项B(" + Str(i) + "):"))

Next i

For i = 1 To n

x(i) = Val(InputBox("请输入初始值X(" + Str(i) + "):"))

Next i

m = 0

f = 0

While f = 0

Max = 0

For i = 1 To n

Sum = 0

For j = 1 To n

Sum = Sum + a(i, j) * x(j)

数值分析上机作业

数值分析上机实验报告 选题:曲线拟合的最小二乘法 指导老师: 专业: 学号: 姓名:

课题八曲线拟合的最小二乘法 一、问题提出 从随机的数据中找出其规律性,给出其近似表达式的问题,在生产实践和科学实验中大量存在,通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中,根据统计数据的含碳量与时间关系,试求含碳量y 与时间t 的拟合曲线。 二、要求 1、用最小二乘法进行曲线拟合; 2、近似解析表达式为()33221t a t a t a t ++=?; 3、打印出拟合函数()t ?,并打印出()j t ?与()j t y 的误差,12,,2,1 =j ; 4、另外选取一个近似表达式,尝试拟合效果的比较; 5、*绘制出曲线拟合图*。 三、目的和意义 1、掌握曲线拟合的最小二乘法; 2、最小二乘法亦可用于解超定线代数方程组; 3、探索拟合函数的选择与拟合精度间的关系。 四、计算公式 对于给定的测量数据(x i ,f i )(i=1,2,…,n ),设函数分布为 ∑==m j j j x a x y 0)()(? 特别的,取)(x j ?为多项式 j j x x =)(? (j=0, 1,…,m )

则根据最小二乘法原理,可以构造泛函 ∑∑==-=n i m j i j j i m x a f a a a H 1 10))((),,,(? 令 0=??k a H (k=0, 1,…,m ) 则可以得到法方程 ???? ??????? ?=????????????????????????),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ????????????????????? 求该解方程组,则可以得到解m a a a ,,,10 ,因此可得到数据的最小二乘解 ∑=≈m j j j x a x f 0)()(? 曲线拟合:实际工作中,变量间未必都有线性关系,如服药后血药浓度与时间的关系;疾病疗效与疗程长短的关系;毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据,并用拟合的曲线方程分析两变量间的关系。 五、结构程序设计 在程序结构方面主要是按照顺序结构进行设计,在进行曲线的拟合时,为了进行比较,在程序设计中,直接调用了最小二乘法的拟合函数polyfit ,并且依次调用了plot 、figure 、hold on 函数进行图象的绘制,最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差,进行拟合效果的比较。

数值计算方法大作业

目录 第一章非线性方程求根 (3) 1.1迭代法 (3) 1.2牛顿法 (4) 1.3弦截法 (5) 1.4二分法 (6) 第二章插值 (7) 2.1线性插值 (7) 2.2二次插值 (8) 2.3拉格朗日插值 (9) 2.4分段线性插值 (10) 2.5分段二次插值 (11) 第三章数值积分 (13) 3.1复化矩形积分法 (13) 3.2复化梯形积分法 (14) 3.3辛普森积分法 (15) 3.4变步长梯形积分法 (16) 第四章线性方程组数值法 (17) 4.1约当消去法 (17) 4.2高斯消去法 (18) 4.3三角分解法 (20)

4.4雅可比迭代法 (21) 4.5高斯—赛德尔迭代法 (23) 第五章常积分方程数值法 (25) 5.1显示欧拉公式法 (25) 5.2欧拉公式预测校正法 (26) 5.3改进欧拉公式法 (27) 5.4四阶龙格—库塔法 (28)

数值计算方法 第一章非线性方程求根 1.1迭代法 程序代码: Private Sub Command1_Click() x0 = Val(InputBox("请输入初始值x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = (Exp(2 * x0) - x0) / 5 If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例:求f(x)=e2x-6x=0在x=0.5附近的根(ep=10-10)

1.2牛顿法 程序代码: Private Sub Command1_Click() b = Val(InputBox("请输入被开方数x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = x0 - (x0 ^ 2 - b) / (2 * b) If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例:求56的值。(ep=10-10)

东南大学数值分析上机题答案

数值分析上机题 第一章 17.(上机题)舍入误差与有效数 设∑=-= N j N j S 2 2 11 ,其精确值为)111-23(21+-N N 。 (1)编制按从大到小的顺序1 -1 ···1-311-21222N S N +++=,计算N S 的通用 程序; (2)编制按从小到大的顺序1 21 ···1)1(111 222-++--+ -=N N S N ,计算N S 的通用程序; (3)按两种顺序分别计算210S ,410S ,610S ,并指出有效位数(编制程序时用单精度); (4)通过本上机题,你明白了什么? 解: 程序: (1)从大到小的顺序计算1 -1 ···1-311-21222N S N +++= : function sn1=fromlarge(n) %从大到小计算sn1 format long ; sn1=single(0); for m=2:1:n sn1=sn1+1/(m^2-1); end end (2)从小到大计算1 21 ···1)1(111 2 22 -++--+-= N N S N function sn2=fromsmall(n) %从小到大计算sn2 format long ; sn2=single(0); for m=n:-1:2 sn2=sn2+1/(m^2-1); end end (3) 总的编程程序为: function p203()

clear all format long; n=input('please enter a number as the n:') sn=1/2*(3/2-1/n-1/(n+1));%精确值为sn fprintf('精确值为%f\n',sn); sn1=fromlarge(n); fprintf('从大到小计算的值为%f\n',sn1); sn2=fromsmall(n); fprintf('从小到大计算的值为%f\n',sn2); function sn1=fromlarge(n) %从大到小计算sn1 format long; sn1=single(0); for m=2:1:n sn1=sn1+1/(m^2-1); end end function sn2=fromsmall(n) %从小到大计算sn2 format long; sn2=single(0); for m=n:-1:2 sn2=sn2+1/(m^2-1); end end end 运行结果:

数值分析上机题目详解

第一章 一、题目 设∑ =-= N N j S 2 j 2 1 1,其精确值为)11 123(21+--N N 。 1) 编制按从大到小的顺序1 1 13112122 2-+??+-+-=N S N ,计算S N 的通用程序。 2) 编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ,计算S N 的通用程序。 3) 按两种顺序分别计算64210,10,10S S S ,并指出有效位数。(编制程序时用单精度) 4) 通过本次上机题,你明白了什么? 二、通用程序 N=input('Please Input an N (N>1):'); AccurateValue=single((0-1/(N+1)-1/N+3/2)/2); Sn1=single(0); for a=2:N; Sn1=Sn1+1/(a^2-1); end Sn2=single(0); for a=2:N; Sn2=Sn2+1/((N-a+2)^2-1); end fprintf('The value of Sn (N=%d)\n',N); fprintf('Accurate Calculation %f\n',AccurateValue); fprintf('Caculate from large to small %f\n',Sn1); fprintf('Caculate from small to large %f\n',Sn2); disp('____________________________________________________')

三、结果 从结果可以看出有效位数是6位。 感想:可以得出,算法对误差的传播有一定的影响,在计算时选一种好的算法可以使结果更为精确。从以上的结果可以看到从大到小的顺序导致大数吃小数的现象,容易产生较大的误差,求和运算从小数到大数所得到的结果才比较准确。

数值分析上机题目

数值分析上机题目 1、 分别用不动点迭代与Newton 法求解方程250x x e -+=的正根与负根。 2、 Use each of the following methods to find a solution in [0.1,1] accurate to within 10^-4 for 4326005502002010x x x x -+--= a. Bisection method b. Newton’s method c. Secant method d. Method of False Position e. Muller’s method 3、 应用Newton 法求f (x )的零点,e=10^-6,这里f (x )=x-sin (x )。 再用求重根的两种方法求f (x )的零点。 4、 应用Newton 法求f (x )的零点,e=10^-6,f(x)=x-sin(x) 再用Steffensen’s method 加速其收敛。 5、 用Neville’s 迭代差值算法,对于函数2 1 (),11125f x x x = -≤≤+进行lagrange 插值。取不同的等分数n=5,10,将区间[-1,1]n 等分,取等距节点。把f(x)和插值多项式的曲线画在同一张图上进行比较。 6、 画狗的轮廓图 7、 Use Romberg integration to compute the following approximations to ? a 、 Determine R1,1,R2,1,R3,1,R4,1and R5,1,and use these approximations to predict the value of the integral. b 、 Determine R2,2 ,R3,3 ,R4,4 ,and R5,5,and modify your prediction. c 、 Determine R6,1 ,R6,2 ,R6,3 ,R6,4 ,R6,5 and R6,6,and modify your prediction.

数值分析上机作业

昆明理工大学工科研究生《数值分析》上机实验 学院:材料科学与工程学院 专业:材料物理与化学 学号:2011230024 姓名: 郑录 任课教师:胡杰

P277-E1 1.已知矩阵A= 10787 7565 86109 75910 ?? ?? ?? ?? ?? ??,B= 23456 44567 03678 00289 00010 ?? ?? ?? ?? ?? ?? ?? ?? ,错误!未找到引用源。 = 11/21/31/41/51/6 1/21/31/41/51/61/7 1/31/41/51/61/71/8 1/41/51/61/71/81/9 1/51/61/71/81/91/10 1/61/71/81/91/101/11?????????????????? (1)用MA TLAB函数“eig”求矩阵全部特征值。 (2)用基本QR算法求全部特征值(可用MA TLAB函数“qr”实现矩阵的QR分解)。解:MA TLAB程序如下: 求矩阵A的特征值: clear; A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; E=eig(A) 输出结果: 求矩阵B的特征值: clear; B=[2 3 4 5 6;4 4 5 6 7;0 3 6 7 8;0 0 2 8 9;0 0 0 1 0]; E=eig(B) 输出结果:

求矩阵错误!未找到引用源。的特征值: clear; 错误!未找到引用源。=[1 1/2 1/3 1/4 1/5 1/6; 1/2 1/3 1/4 1/5 1/6 1/7; 1/3 1/4 1/5 1/6 1/7 1/8; 1/4 1/5 1/6 1/7 1/8 1/9;1/5 1/6 1/7 1/8 1/9 1/10; 1/6 1/7 1/8 1/9 1/10 1/11]; E=eig(错误!未找到引用源。) 输出结果: (2)A= 10 7877565861097 5 9 10 第一步:A0=hess(A);[Q0,R0]=qr(A0);A1=R0*Q0 返回得到: 第二部:[Q1,R1]=qr(A1);A2=R1*Q1

数值计算方法I上机实验考试题

数值计算方法I 上机实验考试题(两题任选一题) 1.小型火箭初始质量为900千克,其中包括600千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉,由此产生30000牛顿的恒定推力.当燃料用尽时引擎关闭。设火箭上升的整个过程中,空气阻力与速度平方成正比,比例系数为0.4(千克/米).重力加速度取9.8米/秒2. A. 建立火箭升空过程的数学模型(微分方程); B. 求引擎关闭瞬间火箭的高度、速度、加速度,及火箭到达最高点的时间和高度. 2.小型火箭初始质量为1200千克,其中包括900千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉,由此产生40000牛顿的恒定推力.当燃料用尽时引擎关闭。设火箭上升的整个过程中,空气阻力与速度平方成正比,比例系数记作k ,火箭升空过程的数学模型为 0)0(,0,01222==≤≤-+?? ? ??-==t dt dx x t t mg T dt dx k dt x d m 其中)(t x 为火箭在时刻t 的高度,m =1200-15t 为火箭在时刻t 的质量,T (=30000牛顿)为推力,g (=9.8米/秒2)为重力加速度, t 1 (=900/15=60秒)为引擎关闭时刻. 今测得一组数据如下(t ~时间(秒),x ~高度(米),v ~速度(米/秒)): 现有两种估计比例系数k 的方法: 1.用每一个数据(t,x,v )计算一个k 的估计值(共11个),再用它们来估计k 。 2.用这组数据拟合一个k . 请你分别用这两种方法给出k 的估计值,对方法进行评价,并且回答,能否认为空气阻力系数k=0.5(说明理由).

数值计算方法作业

数值计算方法作业 姓名:李琦 学号:062410124 求 013=--x x 在x=1.5附近的一个根。 一.牛顿下山法: #include #include float f(float x) /* 定义函数f(x) */ { return x*x*x-x-1; } void main() { float x0,x1=1.5; x0=1; for(;;) { printf (" x0=%f",x0); printf (" x1=%f\n",x1); x0=x1; x1=x0-((x0*x0*x0-x0-1)/(3*x0*x0-1)); if(x0==x1) break; } printf(" x=%f\n",x1); }

二.加权法 #include #include float f(float x) /* 定义函数f(x) */ { return x*x*x-1; } float f1(float x) /* 定义函数f(x)的导数*/ { return 3*x*x; } void main() { float x0,x1=1.5,c; c=f1(x1);x0=1; printf("c=%f\n",c); for(;;) { printf (" x0=%f",x0); printf (" x1=%f\n",x1); x0=x1; x1=(f(x0)-c*x0)/(1-c); if(x0==x1) break; } printf("x=%f\n",x1); }

三.单点弦法: #include #include float f(float x) /* 定义函数f(x) */ { return x*x*x-x-1; } void main() { float x1,x0=1.5,a; a=f(x0); x1=1; for(;;) { printf (" x0=%f",x0); printf (" x1=%f\n",x1); x0=x1; x1=x0-(f(x0)*(x0-1.5)/(f(x0)-a)); if(x0==x1) break; } printf(" x=%f\n",x1); }

《数值计算方法》上机实验报告

《数值计算方法》上机实验报告华北电力大学 实验名称数值il?算方法》上机实验课程名称数值计算方法专业班级:电力实08学生姓名:李超然学号:200801001008 成绩: 指导教师:郝育黔老师实验日期:2010年04月华北电力大学实验报告数值计算方法上机实验报吿一. 各算法的算法原理及计算机程序框图1、牛顿法求解非线性方程 *对于非线性方程,若已知根的一个近似值,将在处展开成一阶 xxfx ()0, fx ()xkk 泰勒公式 "f 0 / 2 八八,fxfxfxxxxx 0 0 0 0 0 kkkk2! 忽略高次项,有 ,fxfxfxxx 0 ()()(),,, kkk 右端是直线方程,用这个直线方程来近似非线性方程。将非线性方程的 **根代入,即fx ()0, X ,* fxfxxx 0 0 0 0, ,, kkk fx 0 fx 0 0,

解出 fX 0 *k XX,, k' fx 0 k 水将右端取为,则是比更接近于的近似值,即xxxxk, Ik, Ik fx ()k 八XX, Ikk* fx()k 这就是牛顿迭代公式。 ,2,计算机程序框图:,见, ,3,输入变量、输出变量说明: X输入变量:迭代初值,迭代精度,迭代最大次数,\0 输出变量:当前迭代次数,当前迭代值xkl ,4,具体算例及求解结果: 2/16 华北电力大学实验报吿 开始 读入 l>k /fx()0?,0 fx 0 Oxx,,01* fx ()0 XX,,,?10 kk, ,1,kN, ?xx, 10 输出迭代输出X输出奇异标志1失败标志

,3,输入变量、输出变量说明: 结束 例:导出计算的牛顿迭代公式,并il ?算。(课本P39例2-16) 115cc (0), 求解结果: 10. 750000 10.723837 10. 723805 10. 723805 2、列主元素消去法求解线性方程组,1,算法原理: 高斯消去法是利用现行方程组初等变换中的一种变换,即用一个不为零的数乘 -个 方程后加只另一个方程,使方程组变成同解的上三角方程组,然后再自下而上 对上三角 3/16 华北电力大学实验报告方程组求解。 列选主元是当高斯消元到第步时,从列的以下(包括)的各元素中选出绝 aakkkkkk 对值最大的,然后通过行交换将其交换到的位置上。交换系数矩阵中的 两行(包括常ekk 数项),只相当于两个方程的位置交换了,因此,列选主元不影响求解的结 ,2,计算机程序框图:,见下页, 输入变量:系数矩阵元素,常向量元素baiji 输出变量:解向量元素bbb,,12n

(完整版)数值计算方法上机实习题答案

1. 设?+=1 05dx x x I n n , (1) 由递推公式n I I n n 1 51+-=-,从0I 的几个近似值出发,计算20I ; 解:易得:0I =ln6-ln5=0.1823, 程序为: I=0.182; for n=1:20 I=(-5)*I+1/n; end I 输出结果为:20I = -3.0666e+010 (2) 粗糙估计20I ,用n I I n n 51 5111+- =--,计算0I ; 因为 0095.05 6 0079.01020 201 020 ≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2 1 20=+= I 程序为:I=0.0087; for n=1:20 I=(-1/5)*I+1/(5*n); end I 0I = 0.0083 (3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。 首先分析两种递推式的误差;设第一递推式中开始时的误差为000I I E '-=,递推过程的舍入误差不计。并记n n n I I E '-=,则有01)5(5E E E n n n -==-=-Λ。因为=20E 20020)5(I E >>-,所此递推式不可靠。而在第二种递推式中n n E E E )5 1(5110-==-=Λ,误差在缩小, 所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中,考虑误差是否得到控制, 即算法是否数值稳定。 2. 求方程0210=-+x e x 的近似根,要求4 1105-+?<-k k x x ,并比较计算量。 (1) 在[0,1]上用二分法; 程序:a=0;b=1.0; while abs(b-a)>5*1e-4 c=(b+a)/2;

东南大学《数值分析》-上机题

数值分析上机题1 设2 21 1N N j S j ==-∑ ,其精确值为1311221N N ??-- ?+?? 。 (1)编制按从大到小的顺序222 111 21311 N S N = +++---,计算N S 的通用程序。 (2)编制按从小到大的顺序22 21111(1)121 N S N N =+++----,计算N S 的通用程序。 (3)按两种顺序分别计算210S ,410S ,610S ,并指出有效位数。(编制程序时用单精度) (4)通过本上机题,你明白了什么? 程序代码(matlab 编程): clc clear a=single(1./([2:10^7].^2-1)); S1(1)=single(0); S1(2)=1/(2^2-1); for N=3:10^2 S1(N)=a(1); for i=2:N-1 S1(N)=S1(N)+a(i); end end S2(1)=single(0); S2(2)=1/(2^2-1); for N=3:10^2 S2(N)=a(N-1); for i=linspace(N-2,1,N-2) S2(N)=S2(N)+a(i); end end S1表示按从大到小的顺序的S N S2表示按从小到大的顺序的S N 计算结果

通过本上机题,看出按两种不同的顺序计算的结果是不相同的,按从大到小的顺序计算的值与精确值有较大的误差,而按从小到大的顺序计算的值与精确值吻合。从大到小的顺序计算得到的结果的有效位数少。计算机在进行数值计算时会出现“大数吃小数”的现象,导致计算结果的精度有所降低,我们在计算机中进行同号数的加法时,采用绝对值较小者先加的算法,其结果的相对误差较小。

东南大学数值分析上机作业汇总

东南大学数值分析上机作业 汇总 -标准化文件发布号:(9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII

数值分析上机报告 院系: 学号: 姓名:

目录 作业1、舍入误差与有效数 (1) 1、函数文件cxdd.m (1) 2、函数文件cddx.m (1) 3、两种方法有效位数对比 (1) 4、心得 (2) 作业2、Newton迭代法 (2) 1、通用程序函数文件 (3) 2、局部收敛性 (4) (1)最大δ值文件 (4) (2)验证局部收敛性 (4) 3、心得 (6) 作业3、列主元素Gauss消去法 (7) 1、列主元Gauss消去法的通用程序 (7) 2、解题中线性方程组 (7) 3、心得 (9) 作业4、三次样条插值函数 (10) 1、第一型三次样条插值函数通用程序: (10) 2、数据输入及计算结果 (12)

作业1、舍入误差与有效数 设∑ =-=N j N j S 2 2 11 ,其精确值为?? ? ??---1112321N N . (1)编制按从小到大的顺序1 1 131121222-? ??+-+-=N S N ,计算N S 的通用程序; (2)编制按从大到小的顺序()1 21 11111222-???+--+-=N N S N ,计算N S 的通用程序; (3)按两种顺序分别计算642101010,,S S S ,并指出有效位数; (4)通过本上机你明白了什么? 程序: 1、函数文件cxdd.m function S=cxdd(N) S=0; i=2.0; while (i<=N) S=S+1.0/(i*i-1); i=i+1; end script 运行结果(省略>>): S=cxdd(80) S= 0.737577 2、函数文件cddx.m function S=cddx (N) S=0; for i=N:-1:2 S=S+1/(i*i-1); end script 运行结果(省略>>): S=cddx(80) S= 0.737577 3、两种方法有效位数对比

(完整版)哈工大-数值分析上机实验报告

实验报告一 题目:非线性方程求解 摘要:非线性方程的解析解通常很难给出,因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。 前言:(目的和意义) 掌握二分法与Newton法的基本原理和应用。 数学原理: 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法:二分法和Newton法。 对于二分法,其数学实质就是说对于给定的待求解的方程f(x),其在[a,b]上连续,f(a)f(b)<0,且f(x)在[a,b]内仅有一个实根x*,取区间中点c,若,则c恰为其根,否则根据f(a)f(c)<0是否成立判断根在区间[a,c]和[c,b]中的哪一个,从而得出新区间,仍称为[a,b]。重复运行计算,直至满足精度为止。这就是二分法的计算思想。

Newton法通常预先要给出一个猜测初值x0,然后根据其迭代公式 产生逼近解x*的迭代数列{x k},这就是Newton法的思想。当x0接近x*时收敛很快,但是当x0选择不好时,可能会发散,因此初值的选取很重要。另外,若将该迭代公式改进为 其中r为要求的方程的根的重数,这就是改进的Newton法,当求解已知重数的方程的根时,在同种条件下其收敛速度要比Newton法快的多。 程序设计: 本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式,如下 function y=f(x); y=-x*x-sin(x); 写成如上形式即可,下面给出主程序。 二分法源程序: clear %%%给定求解区间 b=1.5; a=0;

%%%误差 R=1; k=0;%迭代次数初值 while (R>5e-6) ; c=(a+b)/2; if f12(a)*f12(c)>0; a=c; else b=c; end R=b-a;%求出误差 k=k+1; end x=c%给出解 Newton法及改进的Newton法源程序:clear %%%% 输入函数 f=input('请输入需要求解函数>>','s') %%%求解f(x)的导数 df=diff(f);

数值分析作业

第二章 1. 题目:运用MATLAB编程实现牛顿迭代 2. 实验操作 1、打开MATLAB程序软件。 2、在MATLAB中编辑如下的M程序。 function [p1,err,k,y]=newton(f,df,p0,delta,max) %f 是要求根的方程(f(x)=0); %df 是f(x)的导数; %p0是所给初值,位于x*附近; %delta是给定允许误差; %max是迭代的最大次数; %p1是newton法求得的方程的近似解; %err是p0的误差估计; %k是迭代次数; p0 for k=1:max p1=p0-feval('f',p0)/feval('df',p0); err=abs(p1-p0); p0=p1; k p1 err y=feval('f',p1) if (err> newton('f','df',1.2,10^(-6),20) 3.实验结果

p0 = 1.2000 k =1 p1=1.1030 err=0.0970 y=0.0329 k= 2 p1=1.0524 err=0.0507 y=0.0084 k =3 p1=1.0264 err=0.0260 y=0.0021 k =4 p1=1.0133 err=0.0131 y=5.2963e-004 k =5 p1=1.0066 err=0.0066 y=1.3270e-004 k =6 p1=1.0033 err=0.0033 y=3.3211e-005 k =7 p1=1.0017 err=0.0017 y=8.3074e-006 k =8 p1=1.0008 err=8.3157e-004 y = 2.0774e-006 k =9 p1=1.0004 err=4.1596e-004 y =5.1943e-007 k=10 p1=1.0002 err=2.0802e-004 y= 1.2987e-007 k=11 p1=1.0001 err=1.0402e-004 y =3.2468e-008 k=12 p1=1.0001 err=5.2014e-005 y=8.1170e-009 k=13 p1=1.0000 err=2.6008e-005 y= 2.0293e-009 k=14 p1=1.0000 err=1.3004e-005 y=5.0732e-010 k=15 p1 =1.0000 err=6.5020e-006 y=1.2683e-010 k=16 p1 =1.0000 err=3.2510e-006 y=3.1708e-011 k=17 p1 =1.0000 err=1.6255e-006 y =7.9272e-012 k=18 p1 =1.0000 err =8.1279e-007 y= 1.9820e-012 ans = 1.0000 结果说明:经过18次迭代得到精确解为1,误差为8.1279e-007。

数值计算方法上机实习题

数值计算方法上机实习题 1. 设?+=1 05dx x x I n n , (1) 由递推公式n I I n n 1 51+ -=-,从I 0=0.1824, 0=0.1823I 出发,计算20I ; (2) 20=0I ,20=10000I , 用n I I n n 51 5111+- =--,计算0I ; (3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。 答:第一个算法可得出 e 0=|I 0?I 0 ?| e n =|I n ?I n ?|=5n |e 0| 易知第一个算法每一步计算都把误差放大了5倍,n 次计算后更是放大了5n 倍,可靠性低。 第二个算法可得出 e n =|I n ?I n ?| e 0=(15 )n |e n | 可以看出第二个算法每一步计算就把误差缩小5倍,n 次后缩小了5n 倍,可靠性高。

2. 求方程0210=-+x e x 的近似根,要求41105-+?<-k k x x ,并比较计算量。 (1) 在[0,1]上用二分法; 计算根与步数程序: fplot(@(x) exp(x)+10*x-2,[0,1]); grid on; syms x; f=exp(x)+10*x-2; [root,n]=EFF3(f,0,1); fprintf('root=%6.8f ,n=%d \n',root,n); 计算结果显示: root=0.09057617 ,n=11 (2) 取初值00=x ,并用迭代10 21 x k e x -=+;

(3) 加速迭代的结果; (4) 取初值00 x ,并用牛顿迭代法;

数值计算方法第4次作业

第四章 问题一 一、问题综述 在离地球表面高度为y处的重力加速度如下: 计算高度y=55000m处的重力加速度值。 二、问题分析 以高度y作为自变量,重力加速度的值为因变量。得到以下信息: f(0)=9.8100; f(30000)=9.7487; f(60000)=9.6879; f(90000)=9.6278; f(120000)=9.5682; 本题要求的就是f(55000)的值。 以下将采用课堂中学到的Lagrange插值多项式法、Newton插值多项式法、分段低次插值法和样条插值法求解该问题。 三、问题解决 1. lagrange插值多项式法 对某个多项式函数,已知有给定的k+ 1个取值点: 其中对应着自变量的位置,而对应着函数在这个位置的取值。 假设任意两个不同的x j都互不相同,那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为:

其中每个为拉格朗日基本多项式(或称插值基函数),其表达式为: 拉格朗日基本多项式的特点是在上取值为1,在其它的点上取值为0。 源程序lagrange.m function [c,f]=lagrange(x,y,a) % 输入:x是自变量的矩阵;y是因变量的矩阵;a是要计算的值的自变量; % 输出:c是插值多项式系数矩阵;f是所求自变量对应的因变量; m=length(x); l=zeros(m,m); % l是权矩阵 f=0; for i=1:m v=1; for j=1:m if i~=j v=conv(v,poly(x(j)))/(x(i)-x(j)); % v是l_i(x)的系数矩阵 end end l(i,:)=v; % l矩阵的每一行都是x从高次到低次的系数矩阵 end c=vpa(y*l,10); % 对应阶次的系数相加,乘以y,显示10位有效数字 for k=1:m f=f+c(k)*a^(m-k); end 输入矩阵 x=[0 30000 60000 90000 120000] y=[9.81 9.7487 9.6879 9.6278 9.5682] a=55000 再运行源函数,可得: c = [ -2.057613169e-23, 4.938271605e-18, -3.703703702e-14, -0.000002046111111, 9.81] f = 9.6979851723251649906109417384537

数值分析上机题参考答案.docx

如有帮助欢迎下载支持 数值分析上机题 姓名:陈作添 学号: 040816 习题 1 20.(上机题)舍入误差与有效数 N 1 1 3 1 1 设 S N ,其精确值为 。 2 2 2 N N 1 j 2 j 1 (1)编制按从大到小的顺序 1 1 1 ,计算 S 的通用程序。 S N 1 32 1 N 2 1 N 2 2 (2)编制按从小到大的顺序 1 1 1 ,计算 S 的通用程序。 S N 1 (N 1)2 1 22 1 N N 2 (3)按两种顺序分别计算 S 102 , S 104 , S 106 ,并指出有效位数。 (编制程序时用单精度) (4)通过本上机题,你明白了什么? 按从大到小的顺序计算 S N 的通用程序为: 按从小到大的顺序计算 S N 的通用程序为: #include #include float sum(float N) float sum(float N) { { float j,s,sum=0; float j,s,sum=0; for(j=2;j<=N;j++) for(j=N;j>=2;j--) { { s=1/(j*j-1); s=1/(j*j-1); sum+=s; sum+=s; } } return sum; return sum; } } 从大到小的顺序的值 从小到大的顺序的值 精确值 有效位数 从大到小 从小到大 0.740049 0.74005 0.740049 6 5 S 102 0.749852 0.7499 0.7499 4 4 S 104 0.749852 0.749999 0.749999 3 6 S 106 通过本上机题, 看出按两种不同的顺序计算的结果是不相同的, 按从大到小的顺序计算 的值与精确值有较大的误差, 而按从小到大的顺序计算的值与精确值吻合。 从大到小的顺序 计算得到的结果的有效位数少。 计算机在进行数值计算时会出现“大数吃小数”的现象,导 致计算结果的精度有所降低, 我们在计算机中进行同号数的加法时, 采用绝对值较小者先加 的算法,其结果的相对误差较小。

东南大学-数值分析上机题作业-MATLAB版

2015.1.9 上机作业题报告 JONMMX 2000

1.Chapter 1 1.1题目 设S N =∑1j 2?1 N j=2 ,其精确值为 )1 1 123(21+--N N 。 (1)编制按从大到小的顺序1 1 131121222-+ ??+-+-=N S N ,计算S N 的通用程序。 (2)编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ,计算S N 的通用程序。 (3)按两种顺序分别计算64210,10,10S S S ,并指出有效位数。(编制程序时用单精度) (4)通过本次上机题,你明白了什么? 1.2程序 1.3运行结果

1.4结果分析 按从大到小的顺序,有效位数分别为:6,4,3。 按从小到大的顺序,有效位数分别为:5,6,6。 可以看出,不同的算法造成的误差限是不同的,好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时,误差限随着N 的增大而增大,可见在累加的过程中,误差在放大,造成结果的误差较大。因此,采取从小到大的顺序累加得到的结果更加精确。 2.Chapter 2 2.1题目 (1)给定初值0x 及容许误差ε,编制牛顿法解方程f(x)=0的通用程序。 (2)给定方程03 )(3 =-=x x x f ,易知其有三个根3,0,3321= *=*-=*x x x ○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时,Newton 迭代序列收敛于根x2*。试确定尽可能大的δ。 ○2试取若干初始值,观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。 (3)通过本上机题,你明白了什么? 2.2程序

数值计算大作业

数值计算大作业 题目一、非线性方程求根 1.题目 假设人口随时间和当时人口数目成比例连续增长,在此假设下人口在短期内的增长建立数学模型。 (1)如果令()N t 表示在t 时刻的人口数目,β 表示固定的人口出生率,则人口数目满足微分方程() ()dN t N t dt β=,此方程的解为0()=t N t N e β; (2)如果允许移民移入且速率为恒定的v ,则微分方程变成() ()dN t N t v dt β=+, 此方程的解为 0()=+ (1) t t v N t N e e βββ -; 假设某地区初始有1000000人,在第一年有435000人移入,又假设在第一年年底该地区人口数量1564000人,试通过下面的方程确定人口出生率β,精确到 410-;且通过这个数值来预测第二年年末的人口数,假设移民速度v 保持不变。 435000 1564000=1000000(1) e e βββ + - 2.数学原理 采用牛顿迭代法,牛顿迭代法的数学原理是,对于方程0)(=x f ,如果) (x f 是线性函数,则它的求根是很容易的,牛顿迭代法实质上是一种线性化方法,其基本思想是将非线性方程0)(=x f 逐步归结为某种线性方程来求解。 设已知方程0)(=x f 有近似根k x (假定0)(≠'x f ),将函数)(x f 在点k x 进行泰勒展开,有 . ))(()()(???+-'+≈k k k x x x f x f x f 于是方程0)(=x f 可近似地表示为 ))(()(=-'+k k x x x f x f 这是个线性方程,记其根为1k x +,则1k x +的计算公式为

数值分析上机作业1-1

数值计算方法上机题目1 1、实验1. 病态问题 实验目的: 算法有“优”与“劣”之分,问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化的比较敏感,反之属于好问题。希望读者通过本实验对此有一个初步的体会。 数值分析的大部分研究课题中,如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决,当然一般要付出一些代价(如耗用更多的机器时间、占用更多的存储空间等)。 问题提出: 考虑一个高次的代数多项式 ∏=-= ---=20 1)()20)...(2)(1()(k k x x x x x p (E1-1) 显然该多项式的全部根为l ,2,…,20,共计20个,且每个根都是单重的(也称为简 单的)。现考虑该多项式方程的一个扰动 0)(19 =+x x p ε (E1-2) 其中ε是一个非常小的数。这相当于是对(E1-1)中19 x 的系数作一个小的扰动。我们希望比较(E1-1)和(E1-2)根的差别,从而分析方程(E1-1)的解对扰动的敏感性。 实验内容: 为了实现方便,我们先介绍两个 Matlab 函数:“roots ”和“poly ”,输入函数 u =roots (a ) 其中若变量a 存储1+n 维的向量,则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,...,,+n a a a ,则输出u 的各分量是多项式方程 0...1121=++++-n n n n a x a x a x a 的全部根,而函数 b=poly(v) 的输出b 是一个n +1维变量,它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“Poly ”是两个互逆的运算函数. ve=zeros(1,21); ve(2)=ess; roots(poly(1:20))+ve) 上述简单的Matlab 程序便得到(E1-2)的全部根,程序中的“ess ”即是(E1-2)中的ε。 实验要求: (1)选择充分小的ess ,反复进行上述实验,记录结果的变化并分析它们。如果扰动项的系数ε很小,我们自然感觉(E1-1)和(E1-2)的解应当相差很小。计算中你有什么出乎意料的发现?表明有些解关于如此的扰动敏感性如何? (2)将方程(E1-2)中的扰动项改成18 x ε或其他形式,实验中又有怎样的现象出现?

相关文档
相关文档 最新文档