信计091 龚立丽 200900901004
数值方法计算实习题
要求
对教师党员的评价套管和固井爆破片与爆破装置仓库管理基本要求三甲医院都需要复审吗
:1、用Matlab语言或你熟悉的其他算法语言编写程序,使之尽可能具有通用性;2、根据上机计算实践,对所使用的数值方法的特点、性质、有效性、误差和收敛性等方面进行必要的讨论和分析;3、完成计算后写出实验
报告
软件系统测试报告下载sgs报告如何下载关于路面塌陷情况报告535n,sgs报告怎么下载竣工报告下载
,内容包括:课题名称、解决的问题、采用的数值方法、算法程序、数值结果、对实验结果的讨论和分析等;4、特别说明:严禁抄袭,否则一经发现,所有雷同实验报告最多评为及格。
一、下
表
关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf
给出了飞行中鸭子的上部形状的节点数据,试用三次样条插值函数(自然边界条件)和20次Lagrange插值多项式对数据进行插值。用图示出给定的数据,以及
和
。
0.9
1.3
1.9
2.1
2.6
3.0
3.9
4.4
4.7
5.0
6.0
1.3
1.5
1.85
2.1
2.6
2.7
2.4
2.15
2.05
2.1
2.25
7.0
8.0
9.2
10.5
11.3
11.6
12
12.6
13.0
13.3
2.3
2.25
1.95
1.4
0.9
0.7
0.6
0.5
0.4
0.25
解:>> x=[0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 9.2 10.5 11.3 11.6 12 12.6 13.0 13.3];
>> y=[1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25];
(1)三次样条插值法
在MATLAB中编写m文件
function[f,f0]=scyt(x,y,y2_1,y2_N,x0)
%y2_1和y2_N分别为自然边界条件
%插值点x的坐标:x0
syms t;
f=0.0;f0=0.0;
if(length(x)==length(y))
n=length(x);
else
disp('x和y的维数不相等');
return;
end
for i=1:n
if(x(i)<=x0)&&(x(i+1)>=x0)
index=i;
break;
end
end
A=diag(2*ones(1,n));
A(1,2)=1;
A(n,n-1)=1;
u=zeros(n-2,1);
lamda=zeros(n-1,1);
c=zeros(n,1);
for i=2:n-1
u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));
lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));
c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));
A(i,i+1)=u(i-1);
A(i,i-1)=lamda(i);
end
c(1)=3*(y(2)-y(1))/(x(2)-x(1))-(x(2)-x(1))*y2_1/2;
c(n)=3*(y(n)-y(n-1))/(x(n)-x(n-1))-(x(n)-x(n-1))*y2_N/2;
m=zgf(A,c);
h=x(index+1)-x(index);
f=y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/h;
f0=subs(f,'t',x0);
其中的zgf函数为追赶法,其程序为
function x=zgf(A,b)
n = rank(A);
for(i=1:n)
if(A(i,i)==0)
disp('Error: 对角有元素为0!');
return;
end
end;
d = ones(n,1);
a = ones(n-1,1);
c = ones(n-1);
for(i=1:n-1)
a(i,1)=A(i+1,i);
c(i,1)=A(i,i+1);
d(i,1)=A(i,i);
end
d(n,1) = A(n,n);
for(i=2:n)
d(i,1)=d(i,1) - (a(i-1,1)/d(i-1,1))*c(i-1,1);
b(i,1)=b(i,1) - (a(i-1,1)/d(i-1,1))*b(i-1,1);
end
x(n,1) = b(n,1)/d(n,1);
for(i=(n-1):-1:1)
x(i,1) = (b(i,1)-c(i,1)*x(i+1,1))/d(i,1);
end
在MATLAB中输入指令
>> [f,f0]=scyt(x,y,0,0)
f =
1000/729*(27/5*t-1377/100)*(t-39/10)^2+1000/729*(522/25-24/5*t)*(t-3)^2+100/81*(-9/*t+57/)*(39/10-t)^2-100/81*(-/20+3/2*t)*(t-3)^2
f0 =
2.5851 得三次样条插值函数
S(x)= 1000/729*(27/5*x-1377/100)*(x-39/10)^2+1000/729*(522/25-24/5*x)*(x-3)^2+100/81*(-9/*x+57/)*(39/10-x)^2-100/81*(-/20+3/2*x)*(x-3)^2
>> xi=0.9:0.01:13.3;yi=interp1(x,y,xi,'spline');
>> title('试验一--三次样条插值图示')
(2)用拉格朗日法插值
%定义Lagrange程序
function f=Language(x,y,x0)
syms t;
if(length(x)==length(y))
n=length(x);
else
disp('x和y的维数不相等');
return;
end
f=0.0;
for(i=1:n)
l=y(i);
for(j=1:i-1)
l=l*(t-x(j))/(x(i)-x(j));
end;
for(j=i+1:n)
l=l*(t-x(j))/(x(i)-x(j));
end;
f=f+l;
simplify(f);
if(i==n)
if(nargin==3)
f=subs(f,'t',x0);
else
f=collect(f);
f=vpa(f,6);
end
end
end
>> Language(x,y)
ans =
52462.6*t+189995.*t^3-189851.*t^4+136778.*t^5-11.3161*t^12-.277283e-6*t^18+1.18284*t^13-73866.6*t^6+.111076e-4*t^17-.976904e-1*t^14+.427949e-8*t^19-.307453e-10*t^20+30677.6*t^7+2564.20*t^9-9968.98*t^8+.628590e-2*t^15-525.813*t^10-9652.78-.308159e-3*t^16+86.2514*t^11-128683.*t^2
二、已知Wilson矩阵
,且向量
,则方程组
有准确解
。
⑴用Matlab 内部函数求
,
的所有特征值和
;
⑵令
,解方程组
,并求出向量
和
,从理论结果和实际计算结果两方面分析方程组
解的相对误差
与
的相对误差
的关系;
⑶再改变扰动矩阵
(其元素的绝对值不超过0.005),重复第2问。
解:(1)A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; b=[32;23;33;31];
M=det(A)
M =
1
A的所以特征值:
>> D=eig(A)
D =
0.0102
0.8431
3.8581
30.2887
条件数 >> n=30.2887/0.0102
n =
2.9695e+003
所以A的行列式为1,cond(A)2=2.9695e+003
(2) >> B=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 4.99 9 9.98];
>> b=[32;23;33;31];
>> [rank(B),rank([B,b])]
ans =
4 4
>> x=B\b
x =-5.4761
11.4934
-1.4292
2.4838
求
假设X=
>> x=B\b;x1=[1;1;1;1];X=x-x1
X =
-6.4761
10.4934
-2.4292
1.4838
求
>>norm(X)
ans =
12.6552
12.655就是
。
%求
>> norm(X)/norm(x1)
ans =
6.3276
6.3276 即为
>> norm(a)
ans =
0.2244
>> norm(A)
ans =
30.2887
>> norm(a)/norm(A) %求
ans =
0.0074
所以
=0.0074
inv(A) %求A的逆矩阵,下求
继续阅读