数值分析大作业三四五六七.docx
《数值分析大作业三四五六七.docx》由会员分享,可在线阅读,更多相关《数值分析大作业三四五六七.docx(26页珍藏版)》请在冰点文库上搜索。
![数值分析大作业三四五六七.docx](https://file1.bingdoc.com/fileroot1/2023-6/27/cabf90eb-7a6f-4bbe-8d84-6c6acd4e2371/cabf90eb-7a6f-4bbe-8d84-6c6acd4e23711.gif)
数值分析大作业三四五六七
大作业三
1.给定初值
及容许误差
,编制牛顿法解方程f(x)=0的通用程序.
解:
Matlab程序如下:
函数m文件:
functionFu=fu(x)
Fu=x^3/3-x;
end
函数m文件:
functionFu=dfu(x)
Fu=x^2-1;
end
用Newton法求根的通用程序
clear;
x0=input('请输入初值x0:
');
ep=input('请输入容许误差:
');
flag=1;
whileflag==1
x1=x0-fu(x0)/dfu(x0);
ifabs(x1-x0)flag=0;
end
x0=x1;
end
fprintf('方程的一个近似解为:
%f\n',x0);
寻找最大δ值的程序:
clear
eps=input('请输入搜索精度:
');
ep=input('请输入容许误差:
');
flag=1;
k=0;
x0=0;
whileflag==1
sigma=k*eps;
x0=sigma;
k=k+1;
m=0;
flag1=1;
whileflag1==1&&m<=10^3
x1=x0-fu(x0)/dfu(x0);
ifabs(x1-x0)flag1=0;
end
m=m+1;
x0=x1;
end
ifflag1==1||abs(x0)>=ep
flag=0;
end
end
fprintf('最大的sigma值为:
%f\n',sigma);
2.求下列方程的非零根
解:
Matlab程序为:
(1)主程序
clear
clc
formatlong
x0=765;
N=100;
errorlim=10^(-5);
x=x0-f(x0)/subs(df(),x0);
n=1;
whilenx=x0-f(x0)/subs(df(),x0);
ifabs(x-x0)>errorlim
n=n+1;
else
break;
end
x0=x;
end
disp(['迭代次数:
n=',num2str(n)])
disp(['所求非零根:
正根x1=',num2str(x),'负根x2=',num2str(-x)])
(2)子函数非线性函数f
functiony=f(x)
y=log((513+*x)/*x))-x/(1400*;
end
(3)子函数非线性函数的一阶导数df
functiony=df()
symsx1
y=log((513+*x1)/*x1))-x1/(1400*;
y=diff(y);
end
运行结果如下:
迭代次数:
n=5
所求非零根:
正根x1=负根x2=
大作业四
分析:
(1)输出插值多项式。
(2)在区间[-5,5]内均匀插入99个节点,计算这些节点上函数f(x)的近似值,并在同一张图上画出原函数和插值多项式的图形。
(3)观察龙格现象,计算插值函数在各节点处的误差,并画出误差图。
解:
Matlab程序代码如下:
%此函数实现y=1/(1+4*x^2)的n次Newton插值,n由调用函数时指定
%函数输出为插值结果的系数向量(行向量)和插值多项式
function[ty]=func5(n)
x0=linspace(-5,5,n+1)';
y0=1./(1.+4.*x0.^2);
b=zeros(1,n+1);
fori=1:
n+1
s=0;
forj=1:
i
t=1;
fork=1:
i
ifk~=j
t=(x0(j)-x0(k))*t;
end;
end;
s=s+y0(j)/t;
end;
b(i)=s;
end;
t=linspace(0,0,n+1);
fori=1:
n
s=linspace(0,0,n+1);
s(n+1-i:
n+1)=b(i+1).*poly(x0(1:
i));
t=t+s;
end;
t(n+1)=t(n+1)+b
(1);
y=poly2sym(t);
10次插值运行结果:
[bY]=func5(10)
b=
Columns1through4
Columns5through8
Columns9through11
Y=
b为插值多项式系数向量,Y为插值多项式。
插值近似值:
x1=linspace(-5,5,101);
x=x1(2:
100);
y=polyval(b,x)
y=
Columns1through12
Columns13through24
Columns25through36
Columns37through48
Columns49through60
Columns61through72
Columns73through84
Columns85through96
Columns97through99
绘制原函数和拟合多项式的图形代码:
plot(x,1./(1+4.*x.^2))
holdall
plot(x,y,'r')
xlabel('X')
ylabel('Y')
title('Runge现象')
gtext('原函数')
gtext('十次牛顿插值多项式')
绘制结果:
误差计数并绘制误差图:
holdoff
ey=1./(1+4.*x.^2)-y
ey=
Columns1through12
Columns13through24
Columns25through36
Columns37through48
Columns49through60
0
Columns61through72
Columns73through84
Columns85through96
Columns97through99
plot(x,ey)
xlabel('X')
ylabel('ey')
title('Runge现象误差图')
输出结果为:
大作业五
解:
Matlab程序为:
x=[-520,-280,,-78,,,0,,,78,,280,520]';
y=[0,-30,-36,-35,,,0,,,35,36,30,0]';
n=13;
%求解M
fori=1:
1:
n-1
h(i)=x(i+1)-x(i);
end
fori=2:
1:
n-1
a(i)=h(i-1)/(h(i-1)+h(i));
b(i)=1-a(i);
c(i)=6*((y(i+1)-y(i))/h(i)-(y(i)-y(i-1))/h(i-1))/(h(i-1)+h(i));
end
a(n)=h(n-1)/(h
(1)+h(n-1));
b(n)=h
(1)/(h
(1)+h(n-1));
c(n)=6/(h
(1)+h(n-1))*((y
(2)-y
(1))/h
(1)-(y(n)-y(n-1))/h(n-1));
A(1,1)=2;
A(1,2)=b
(2);
A(1,n-1)=a
(2);
A(n-1,n-2)=a(n);
A(n-1,n-1)=2;
A(n-1,1)=b(n);
fori=2:
1:
n-2
A(i,i)=2;
A(i,i+1)=b(i+1);
A(i,i-1)=a(i+1);
end
C=c(2:
n);
C=C';
m=A\C;
M
(1)=m(n-1);
M(2:
n)=m;
xx=-520:
:
520;
fori=1:
51
forj=1:
1:
n-1
ifx(j)<=xx(i)&&xx(i)break;
end
end
yy(i)=M(j+1)*(xx(i)-x(j))^3/(6*h(j))-M(j)*(xx(i)-x(j+1))^3/(6*h(j))+(y(j+1)-M(j+1)*h(j)^2/6)*(xx(i)-x(j))/h(j)-(y(j)-M(j)*h(j)^2/6)*(xx(i)-x(j+1))/h(j);
end;
fori=52:
101
yy(i)=-yy(102-i);
end;
fori=1:
50
xx(i)=-xx(i);
end
plot(xx,yy);
holdon;
fori=1:
1:
n/2
x(i)=-x(i);
end
plot(x,y,'bd');
title('机翼外形曲线');
输出结果:
运行文件,得到
2.
(1)编制求第一型3次样条插值函数的通用程序;
(2)已知汽车门曲线型值点的数据如下:
解:
(1)Matlab编制求第一型3次样条插值函数的通用程序:
function[Sx]=Threch(X,Y,dy0,dyn)
%X为输入变量x的数值
%Y为函数值y的数值
%dy0为左端一阶导数值
%dyn为右端一阶导数值
%Sx为输出的函数表达式
n=length(X)-1;
d=zeros(n+1,1);
h=zeros(1,n-1);
f1=zeros(1,n-1);
f2=zeros(1,n-2);
fori=1:
n%求函数的一阶差商
h(i)=X(i+1)-X(i);
f1(i)=(Y(i+1)-Y(i))/h(i);
end
fori=2:
n%求函数的二阶差商
f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));
d(i)=6*f2(i);
end
d
(1)=6*(f1
(1)-dy0)/h
(1);
d(n+1)=6*(dyn-f1(n-1))/h(n-1);
%赋初值
A=zeros(n+1,n+1);
B=zeros(1,n-1);
C=zeros(1,n-1);
fori=1:
n-1
B(i)=h(i)/(h(i)+h(i+1));
C(i)=1-B(i);
end
A(1,2)=1;
A(n+1,n)=1;
fori=1:
n+1
A(i,i)=2;
end
fori=2:
n
A(i,i-1)=B(i-1);
A(i,i+1)=C(i-1);
end
M=A\d;
symsx;
fori=1:
n
Sx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))+M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);
digits(4);
Sx(i)=vpa(Sx(i));
end
fori=1:
n
disp('S(x)=');
fprintf('%s(%d,%d)\n',char(Sx(i)),X(i),X(i+1));
end
S=zeros(1,n);
fori=1:
n
x=X(i)+;
S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))+M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3;
end
disp('S(i+');
disp('iX(i+S(i+');
fori=1:
n
fprintf('%d%.4f%.4f\n',i,X(i)+,S(i));
end
end
输出结果:
>>X=[012345678910];
>>Y=[];
>>Threch(X,Y,,
S(x)=
*x-*x^2-*x^3+(0,1)
S(x)=
*x-*x^2-*x^3+(1,2)
S(x)=
*x-*x^2-*x^3+(2,3)
S(x)=
*x^2-*x-*x^3+(3,4)
S(x)=
*x-*x^2+*x^3-(4,5)
S(x)=
*x^2-*x-*x^3+(5,6)
S(x)=
*x-*x^2+*x^3-(6,7)
S(x)=
*x^2-*x-*x^3+(7,8)
S(x)=
*x-*x^2+*x^3-(8,9)
S(x)=
*x-*x^2+*x^3-(9,10)
S(i+
iX(i+S(i+
1
2
3
4
5
6
7
8
9
10
ans=
[-*x^3-*x^2+*x+,-*x^3-*x^2+*x+,-*x^3-*x^2+*x+,-*x^3+*x^2-*x+,*x^3-*x^2+*x-,-*x^3+*x^2-*x+,*x^3-*x^2+*x-,-*x^3+*x^2-*x+,*x^3-*x^2+*x-,*x^3-*x^2+*x-]
大作业六
1、炼钢厂出钢时所用的圣刚睡的钢包,在使用过程中由于钢液及炉渣对包衬耐火材料的侵蚀,使其容积不断增大,经试验,钢包的容积与相应的使用次数的数据如下:
(使用次数x,容积y)
x
y
x
y
2
9
3
10
5
14
6
16
7
17
9
19
10
20
选用双曲线
对使用最小二乘法数据进行拟合。
解:
Matlab程序如下:
functiona=nihehanshu()
x0=[2356791011121416171920];
y0=[];
A=zeros(2,2);
B=zeros(2,1);
a=zeros(2,1);
x=1./x0;
y=1./y0;
A(1,1)=14;
A(1,2)=sum(x);
A(2,1)=A(1,2);
A(2,2)=sum(x.^2);
B
(1)=sum(y);
B
(2)=sum(x.*y);
a=A\B;
y=1./(a
(1)+a
(2)*1./x0);
subplot(1,2,2);
plot(x0,y0-y,'bd-');
title('拟合曲线误差');
subplot(1,2,1);
plot(x0,y0,'go');
holdon;
x=2:
:
20;
y=1./(a
(1)+a
(2)*1./x);
plot(x,y,'r*-');
legend('散点','拟合曲线图1/y=a
(1)+a
(2)*1/x');
title('最小二乘法拟合曲线');
试验所求的系数为:
nihehanshu
ans=
则拟合曲线为
拟合曲线图、散点图、误差图如下:
2、下面给出的是乌鲁木齐最近1个月早晨7:
00左右(新疆时间)的天气预报所得到的温度,按照数据找出任意次曲线拟合方程和它的图像。
用MATLAB编程对上述数据进行最小二乘拟合。
2008年10月--11月26日
天数
1
2
3
4
5
6
7
8
9
10
温度
9
10
11
12
13
14
13
12
11
9
天数
11
12
13
14
15
16
17
18
19
20
温度
10
11
12
13
14
12
11
10
9
8
天数
21
22
23
24
25
26
27
28
29
30
温度
7
8
9
11
9
7
6
5
3
1
解:
Matlab的程序如下:
x=[1:
1:
30];
y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1];
a1=polyfit(x,y,3)%三次多项式拟合%
a2=polyfit(x,y,9)%九次多项式拟合%
a3=polyfit(x,y,15)%十五次多项式拟合%
b1=polyval(a1,x)
b2=polyval(a2,x)
b3=polyval(a3,x)
r1=sum((y-b1).^2)%三次多项式误差平方和%
r2=sum((y-b2).^2)%九次次多项式误差平方和%
r3=sum((y-b3).^2)%十五次多项式误差平方和%
plot(x,y,'*')%用*画出x,y图像%
holdon
plot(x,b1,'r')%用红色线画出x,b1图像%
holdon
plot(x,b2,'g')%用绿色线画出x,b2图像%
holdon
plot(x,b3,'b:
o')%用蓝色o线画出x,b3图像%
试验结果为:
a1=
Columns1through2
Columns3through4
a2=
Columns1through2
Columns3through4
Columns5through6
Columns7through8
Columns9through10
a3=
Columns1through2
Columns3through4
Columns5through6
Columns7through8
Columns9through10
Columns11through12
Columns13through14
Columns15through16
b1=
Columns1through2
Columns3through4
Columns5through6
Columns7through8
Columns9through10
Columns11through12
Columns13through14
Columns15through16
Columns17through18
Columns19through20
Columns21through22
Columns23through24
Columns25through26
Columns27through28
Columns29through30
b2=
Columns1through2
Columns3through4
Columns5through6
Columns7through8
Columns9through10
Columns11through12
Columns13through14
Columns15through16
Columns17through18
Columns19through20
Columns21through22
Columns23through24
Columns25through26
Columns27through28
Columns29through30
b3=
Columns1through2
Columns3through4
Columns5through6
Columns7through8
Columns9through10
Columns11through12
Columns13through14
Columns15through16
Columns17through18
Columns19through20
Columns21through22
Columns23through24
Columns25through26
Columns27through28
Columns29through30
r1=
r2=
r3=
其中的最后图像为:
大作业七
用Romberg求积法计算积分
的近似值,要求误差不超过
。
解:
Matlab程序为:
%被积函数m文件:
functionFx=fx(x)
Fx=1/(1+100*x*x);
end
%Romberg求积法计算积分的通用程序
functionRomberg()
clear;
a=input('请输入积分下限:
a=');
b=input('请输入积分上限:
b=');
eps=input('请输入允许精度:
eps=');
%========计算Tn========%
functionTn=T(n)
Tn=0;
h=(b-a)/n;
x=zeros(1,n+1);
fork=1:
n+1
x(k)=a+(k-1)*h;
end
forj=1:
n
Tn=Tn+h*(fx(x(j))+fx(x(j+1)))/2;
end
end
%========计算Sn========%
functionSn=S(n)
Sn=4/3*T(2*n)-1/3*T(n);
end
%========计算Cn========%
functionCn=C(n)
Cn=16/15*S(2*n)-1/15*S(n);
end
%========计算Rn========%
functionRn=R(n)
Rn=64/63*C(2*n)-1/63*C(n);
end
%========计算满足允许精度的Rn,并打印输出========%
i=1;
flag=1;
whileflag==1
ifabs(R(2^i)-R(2^(i-1)))/255flag=0;
end
i=i+1;
end
fprintf('该积分的值为:
%f\n',R(2^(i-1)));
End
运行结果为:
>>Romberg
请输入积分下限:
a=-1
请输入积分上限:
b=1
请输入允许精度:
eps=*10^-7
该积分的值为:
大作业八
设常微分方程初值问题
其精确解为
选取步长h使经典四阶RK法稳定,求解微分方程