《应用计算方法教程》大作业.doc

《应用计算方法教程》作业

姓名： * *
学号： S2017****
专业： * * * *
学院： 机械工程学院
联系方式： ***********
任课教师： 丁军

作业一.用任意方法求解七阶方程-5+22x-12x2+x3+13x4+10x5-7x6-4x7
牛顿迭代法:
y=@(x)4*x^7+7*x^6-10*x^5-13*x^4-x^3+12*x^2-22*x+5;
s=@(m)28*m^6+42*m^5-50*m^4-52*m^3-3*m^2+24*m-22;
p0=0;N=1000;Tol=0.001;
for k=1:N
p1=p0-y(p0)/s(p0);
if abs(p1-p0)>Tol
p0=p1;
else break
end
end
disp(p1);disp(k)

作业二.实验4-2
合理利用LU分解法（使用列选主元可恰当加分）或追赶法求解下列方程组：

程序：
A=[1.3,-1,-1,2.4,-3.4;2.4,-1,-1,1.4,3.7;2,1,-2,3.6,6.8;2.5,-1,4,3,6.6;1.5,-1,-1,5.3,2.8];

b=[4.2;6.3;5.5;3.6;6.2];
L=eye(5);U(1,:)=A(1,:);
y=[0;0;0];x=[0;0;0];
for k=2:5
if U(k-1,k-1)==0
{
disp('分解失败');
Return
}
end
L(k:5,k-1)=A(k:5,k-1)/U(k-1,k-1)
U(k,k:5)=A(k,k:5)-L(k,1:k-1)*U(1:k-1,k:5)
if k<5
A(k+1:5,k)=A(k+1:5,k)-L(k+1:5,1:k-1)*U(1:k-1,k);
end
end
L1=inv(L);y=L1*b;
U1=inv(U);x=U1*y;
disp(L);disp(U);disp(x);

结果：

作业三.实验4-3

1
1/2
1/3
1/4
1/5
1/6
1/7
1/8
1/9
1/10
1/2
1/3
1/4
1/5
1/6
1/7
1/8
1/9
1/10
1/11
1/3
1/4
1/5
1/6
1/7
1/8
1/9
1/10
1/11
1/12
1/4
1/5
1/6
1/7
1/8
1/9
1/10
1/11
1/12
1/13
1/5
1/6
1/7
1/8
1/9
1/10
1/11
1/12
1/13
1/14
1/6
1/7
1/8
1/9
1/10
1/11
1/12
1/13
1/14
1/15
1/7
1/8
1/9
1/10
1/11
1/12
1/13
1/14
1/15
1/16
1/8
1/9
1/10
1/11
1/12
1/13
1/14
1/15
1/16
1/17
1/9
1/10
1/11
1/12
1/13
1/14
1/15
1/16
1/17
1/18

1/10
1/11
1/12
1/13
1/14
1/15
1/16
1/17
1/18
1/19
观察结论：矩阵的条件数越来越大，而且远远大于1，这是一个病态问题

1） LU分解法
算法实现：
% n=5:10
clear;clc;
n=10
for i=1:n
for j=1:n
H(i,j)=1/(i+j-1);
end
end
disp(H)
cond(H,1)
format short;
n=size(H);
L=eye(n);
U(1,:) = H(1,:);
for k = 2:n
if U(k-1,k-1)==0
disp('分解失败');return
end
L(k:n,k-1)=H(k:n,k-1)/U(k-1,k-1);
U(k,k:n)=H(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);
if k<n
H(k+1:n,k)=H(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k);
end
end
B([1:10],1)=1;
Y=L\B;
L
U
X=U\Y
运行结果：
L=
1.0000
0
0
0
0
0
0
0
0
0
0.5000
1.0000
0
0
0
0
0
0
0
0
0.3333
1.0000
1.0000
0
0
0
0
0
0
0
0.2500
0.9000
1.5000
1.0000
0
0
0
0
0
0
0.2000
0.8000
1.7143
2.0000
1.0000
0
0
0
0
0
0.1667
0.7143
1.7857
2.7778
2.5000
1.0000
0
0
0
0
0.1429
0.6429
1.7857
3.3333
4.0909
3.0000
1.0000
0
0
0
0.1250
0.5833
1.7500
3.7121
5.5682
5.6538
3.5000
1.0000
0
0

0.1111
0.5333
1.6970
3.9596
6.8531
8.6154
7.4667
4.0000
1.0000
0
0.1000
0.4909
1.6364
4.1119
7.9301
11.6308
12.6000
9.5294
4.5000
1.0000
U=
1.0000
0.5000
0.3333
0.2500
0.2000
0.1667
0.1429
0.1250
0.1111
0.1000
0
0.0833
0.0833
0.0750
0.0667
0.0595
0.0536
0.0486
0.0444
0.0409
0
0
0.0056
0.0083
0.0095
0.0099
0.0099
0.0097
0.0094
0.0091
0
0
0
0.0004
0.0007
0.0010
0.0012
0.0013
0.0014
0.0015
0
0
0
0
0.0000
0.0001
0.0001
0.0001
0.0002
0.0002
0
0
0
0
0
0.0000
0.0000
0.0000
0.0000
0.0000
0
0
0
0
0
0
0.0000
0.0000
0.0000
0.0000
0
0
0
0
0
0
0
0.0000
0.0000
0.0000
0
0
0
0
0
0
0
0
0.0000
0.0000
0
0
0
0
0
0
0
0
0
0.0000
X=1.0e+06 *
-0.0000
0.0010
-0.0238
0.2402
-1.2611
3.7834
-6.7261
7.0007
-3.9379
0.9237
2） LU分解迭代求精法
算法实现：
% n=5:10
clear;clc;
n=10;
for i=1:n
for j=1:n
H(i,j)=1/(i+j-1);
end
end

format short;
n=size(H);
L=eye(n);
U(1,:) = H(1,:);
tol=10^(-3);
for k = 2:n
if U(k-1,k-1)==0
disp('分解失败');return
end
L(k:n,k-1)=H(k:n,k-1)/U(k-1,k-1);
U(k,k:n)=H(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);
if k<n
H(k+1:n,k)=H(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k);
end
end
B([1:10],1)=1;
Y=L\B;
X=U\Y;
r=B-H*X;
Y=L\r;
w=U\Y;
format long
while norm(w)<tol
r=B-H*X;
Y=L\r;
w=U\Y;
X=X+w;
end
X
运算结果：
X =1.0e+06 *
-0.000009998346355
0.000989857965323
-0.023756990291064
0.240212759410420
-1.261130560338435
3.783425339621443
-6.726139778865695
7.000720697677755
-3.937927027136435
0.923715699219681

结果分析：对于同一复杂程度的方程组，LU分解迭代求精法所得的解较LU分解法所得到的解的精确度更高。

作业四.实验5-1
分别用Jacobi、Seidel、Sor（=0.8,1.1,1.2,1.3,1.4,1.5）迭代求解下面的方程组，并做结果分析。
初值，精度要求：<。

Jacobi:
A=[12.3,-2,-1,3.4,-3.7;1.4,9,-3,2.4,2.7;2.1,1,8,2.6,5.8;3.5,-2.1,1,13,4.6;2.5,-1,-2,5.3,14.8];
b=[4.8;2.3;2.5;3.6;2.2];
X0=[0;0;0;0;0];
X=X0;K=1;
while K<=30
for i=1:5
X(i)=(b(i)-A(i,:)*X0)/A(i,i)+X0(i);
end
if norm(X-X0)/X<0.00001
disp(X);
disp(K);
return;
end
K=K+1;
X0=X;
end
结果：

Seidel:
A=[12.3,-2,-1,3.4,-3.7;1.4,9,-3,2.4,2.7;2.1,1,8,2.6,5.8;3.5,-2.1,1,13,4.6;2.5,-1,-2,5.3,14.8];
b=[4.8;2.3;2.5;3.6;2.2];
X0=[0;0;0;0;0];
X=X0;K=1;
while K<=20
for i=1:5
X(i)=(b(i)-A(i,:)*X)/A(i,i)+X(i);
end
if norm(X-X0)/X<0.00001
disp(X);
disp(K);
return;
end
K=K+1;
X0=X;
end
disp('发散');
结果：

SOR:
A=[12.3,-2,-1,3.4,-3.7;1.4,9,-3,2.4,2.7;2.1,1,8,2.6,5.8;3.5,-2.1,1,13,4.6;2.5,-1,-2,5.3,14.8];
b=[4.8;2.3;2.5;3.6;2.2];
X0=[0;0;0;0;0];
X=X0;K=1;
w=[0.8,1.1,1.2,1.3,1.4,1.5];
for j=1:6
m=w(j);
while K<=20
for i=1:5
X(i)=m*(b(i)-A(i,:)*X)/A(i,i)+X(i);
end
if norm(X-X0,inf)/norm(X,inf)<0.00001
disp(X);

disp(K);
break;
end
K=K+1;
X0=X;
end
end
结果：
w=0.8 w=1.1

w=1.2 w=1.2

作业五.实验7-2
试验一：
1. 算法实现
% 原函数图像

figure(1)
x=-1:0.01:1;
y=1./(1+25.*x.^2);
plot(x,y,'k')
grid on;hold on;
%Lagrange插值
yy=0;n=5;x=-1:2/n:1;
y=1./(1+25.*x.^2);
xx=-1:0.01:1;
for k=1:n+1
t=1;
for i=1:n+1
if i~=k;
t=t.*(xx-x(i))/(x(k)-x(i));
end
end
yy=yy+t*y(k);
end
plot(xx,yy,'g-')
hold on
n=10;yy=0;x=-1:2/n:1;
y=1./(1+25.*x.^2);xx=-1:0.01:1;
for k=1:n+1
t=1;
for i=1:n+1
if i~=k;
t=t.*(xx-x(i))/(x(k)-x(i));
end
end
yy=yy+t*y(k);
end
plot(xx,yy,'r:*');
legend('y=1./(1+25.*x.^2)','n=5','n=10');
title('Lagrange')
2. 运行结果

试验二：
1. 算法实现
% 原函数图像
figure(1)
x=-1:0.01:1;
y=1./(1+25.*x.^2);
plot(x,y,'k')
grid on;hold on;
% chepyshev结点
n=10;yy=0;kk=0:10;
x=-cos((kk*pi)/n);
y=1./(1+25.*x.^2);
xx=-1:0.01:1;
for k=1:n+1
t=1;
for i=1:n+1
if i~=k;
t=t.*(xx-x(i))/(x(k)-x(i));
end
end
yy=yy+t*y(k);
end
plot(xx,yy,'r:*');
legend('y=1./(1+25.*x.^2)','n=5','n=10');
title('Lagrange and Chebyshev')

运行结果
作业六.实验8-1

算法实现：
clear all;clc;
format short
x=[0.4 0.55 0.65 0.80 0.90 1.05];
y=[0.41075 0.57815 0.69675 0.88811 1.02652 1.25386];
A=zeros(4,4);
for m=1:4
for n=1:4
for i=1:6
A(m,n)=A(m,n)+x(i)^(m+n-2);
end
end
end
d=zeros(4,1);
for n=1:4
for i=1:6
d(n,1)=d(n,1)+y(i)*x(i)^(n-1);
end
end
a=A\d;
R=0;
for i=1:6
R=R+(y(i)-(a(1)+a(2)*x(i)+a(3)*x(i)^2+a(4)*x(i)^3))^2;
end
disp('系数a=');disp(a);disp('平方误差R=');disp(R);

运行结果：
系数a
平方误差R
-0.0056
5.0057e-09

1.0331

-0.0660
0.2137

2）

运用命令行窗口：
Step1：
>> syms x;
>> p0=1;
>> p1=x-quad('[x]',0.40,1.05)/(1.05-0.40)*p0
输出：p1 =x - 29/40
Step2：
>>p2=x^2-quad('[x.^2]',0.40,1.05)/(1.05-0.40)*p0-quad('[x.^3-29*x.^2/40]',0.40,1.05)/quad('[(x-29/40).^2]',0.40,1.05)*p1
输出：p2 =x^2 - (29*x)/20 + 1177/2400

Step3：
>> p3=x^3-quad('[x.^3]',0.40,1.05)/(1.05-0.40)*p0-quad('[x.^4-29*x.^3/40]',0.40,1.05)/quad('[(x-29/40).^2]',0.40,1.05)*p1-quad ('[x.^5 - 29*x.^4/20 + 1177/2400*x.^3]',0.40,1.05)/ quad('[(x.^2 - (29*x)/20 + 1177/2400).^2]',0.40,1.05)*p2
输出：p3 =x^3 - (87*x^2)/40 + (3027*x)/2000 - 53621/160000
Step4：
>> x=[0.40,0.55,0.65,0.80,0.90,1.05];y=[0.41075;0.57815;0.69675;0.88811;1.02652;1.25386];
A(1:6,1)=ones(size(x));
A(1:6,2)=x - 29/40 ;
A(1:6,3)=x.^2 - (29*x)/20 + 1177/2400 ;
A(1:6,4)=x.^3 - (87*x.^2)/40 + (3027*x)/2000 - 53621/160000; C=A\y
输出：
C =
0.8042
1.2881
0.3989
0.2137

Step5：
>> s=0;k=1;
while k<7
s=s+(0.2137* (x(k)^3 - (87*x(k)^2)/40 + (3027*x(k))/2000 - 53621/160000) -0.3989*(x(k)^2 - (29*x(k))/20 + 1177/2400 )+1.2881* (x(k) - 29/40)+0.8042-y(k))^2 ;k=k+1;
end
>> s
输出：
s =
0.0075

结果分析：比较(1)和（2）可以发现正交之后得到的最小二乘法的误差是直接多项式拟合的万分之一，通过正交的pj(x)的函数系解决最小二乘法的病态结果问题，提供了求最小二乘法拟合多项式的快捷方法。

作业七.计算
（1）编写复化梯形公式和复化Simpson公式通用子程序，，分别采用4,8,16,32,64等分区间计算。
（2）使用Romberg求积公式。
(1) 复化梯形公式
程序：
a=0;b=1;s=0;
f=@(x)(log(1+x)/x);
fa=1;fb=log(2);
n=2;
for i=1:5
n=n*2;
h=(b-a)/(n);
f1=0;
for i=1:n-1
f1=f1+f(a+h);
end
s=h*(fa+fb+2*f1)/2;
disp(s)
end
结果：

复化simposon公式
程序：
a=0;b=1;s=0;
f=@(x)(log(1+x)/x);
fa=1;fb=log(2);
n=1;
for j=1:5
n=n*2;
h=(b-a)/(2*n);
f1=0;
for i=0:n-1
f1=f1+f(a+(2*i+1)*h);
end
f2=0;
for i=1:n-1
f2=f2+f(a+2*i*h);
end
s=h*(fa+fb+4*f1+2*f2)/3;
disp(s)
end
结果：

(2)程序：
function y=romberg(f,n,a,b)
f=@(x)(log(1+x)/x);
a=0;b=1;n=4;
z=zeros(n,n);
h=b-a;
z(1,1)=(h/2)*(f(a)+f(b));
f1=0;
for
i=2:n
for
k=1:2^(i-2)
f1=f1+f(a+(k-0.5)*h);
end
z(i,1)=0.5*z(i-1,1)+0.5*h*f1;
h=h/2;
f1=0;
for
j=2:i
z(i,j)=z(i,j-1)+(z(i,j-1)-z(i-1,j-1))/(4^(j-1)-1);
end
end
z,n
fun=inline('ln(1+x)./x');
romberg(fun,8,0,1)

作业八.实验10-1
利用Euler法（h=0.025）, 改进Euler法（h=0.05），4阶Rung-Kutta方法（h=0.1）求解初值问题。
（1）Euler法

程序：
f=@(x,y)(-(y+1)*(y+3));
x0=0;
y0=-2;
n=8;
h=0.025;
x=x0;y=y0;
for i=1:n
y=y+h*f(x,y);
x=x+h;
disp(x)
disp(y)
end
结果：

（2）改进Euler法
程序：
f=@(x,y)(-(y+1)*(y+3));
x0=0;
y0=-2;
n=4;
h=0.05;
x=x0;y=y0;

for i=1:n
y1=y+h*f(x,y);
y2=y+h*f(x+h,y1);
y=(y1+y2)/2;
x=x+h;
disp(x)
disp(y)
end
结果：

（3）4阶Rung-Kutta方法
程序：
f=@(x,y)(-(y+1)*(y+3));
x0=0;
y0=-2;
n=3;
h=0.1;
x=x0;y=y0;
for i=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;
disp(x)
disp(y)

end
结果：

（注：专业文档是经验性极强的领域，无法思考和涵盖全面，素材和资料部分来自网络，供参考。可复制、编制，期待你的好评与关注）