数值计算方法实验报告.docx
- 文档编号:13639867
- 上传时间:2023-06-15
- 格式:DOCX
- 页数:13
- 大小:30.22KB
数值计算方法实验报告.docx
《数值计算方法实验报告.docx》由会员分享,可在线阅读,更多相关《数值计算方法实验报告.docx(13页珍藏版)》请在冰点文库上搜索。
数值计算方法实验报告
数值计算方法实验报告
叶耀北师珠2016.01.10
&1:
实验一:
基本变量与运算:
for循环:
s=0
fori=1:
n(n为常数)
s=s+i;
end
s
函数及函数实验方法:
分段函数的例子:
function[y]=func(x)
ifx>0
y=sin(x);
else
y=1-cos(x);
end
基本绘图命令实验:
下面看如何画出这个分段函数的图形:
x=-pi:
0.1:
pi;
y=[];
fori=1:
length(x)
y(i)=func(x(i));
end
plot(x,y)
画图的程序:
x=-5:
0.1:
5;
y=[];
fori=1:
length(x)
y(i)=func2(x(i));
end
plot(x,y)
gridon
&2:
lagrange插值多项式实现:
完善的拉格朗日插值基函数的程序:
function[f]=language(x,y,x0)
%theLanguageploynomialfunction
%thevectorxandyareknown
symst;
if(length(x)==length(y))
n=length(x);
else
disp('xandyhavediffierentdimensions!
')
return;
end
f=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
end
&3:
Newton插值多项式的程序:
function[y,R,A,C,L]=newdscg(X,Y,x,M)
n=length(X);m=length(x);
fort=1:
m
z=x(t);A=zeros(n,n);A(:
1)=Y';
s=0.0;p=1.0;q1=1.0;c1=1.0;
forj=2:
n
fori=j:
n
A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1));
end
q1=abs(q1*(z-X(j-1)));c1=c1*j;
end
C=A(n,n);q1=abs(q1*(z-X(n)));
fork=(n-1):
-1:
1
C=conv(C,poly(X(k)));
d=length(C);C(d)=C(d)+A(k,k);
end
y(k)=polyval(C,z);
end
R=M*q1/c1;L(k,:
)=poly2sym(C)
&4:
龙贝格积分法的上机实现:
formatlong
f=input('请输入原函数f=','s');
a=input('积分下限a=');
b=input('积分上限b=');
eps1=input('精度eps1=');
T
(1)=double((b-a)/2*(limit(sym(f),findsym(sym(f)),a)+limit(sym(f),findsym(sym(f)),b)));
fork=2:
4
sum1=0;
fori=1:
2^(k-2)
sum1=sum1+subs(sym(f),findsym(sym(f)),(a+(2*i-1)*(b-a)/2^(k-1)));
end
T(k)=1/2*T(k-1)+(b-a)/(2^(k-1))*sum1;
end
fork=1:
3
S(k)=T(k+1)+1/(4-1)*(T(k+1)-T(k));
end
fork=1:
2
c(k)=S(k+1)+1/(4^2-1)*(S(k+1)-S(k));
end
R
(1)=C
(2)+1/(4^3-1)*(C
(2)-C
(1));
k=3;
while1
T
(1)=T
(2);
T
(2)=T(3);
T(3)=T(4);
sum2=0;
fori=1:
2^k
sum2=sum2+subs(sym(f),findsym(sym(f)),(a+(2*i-1)*(b-a)/2^(k+1)));
end
T(4)=1/2*T(4)+(b-a)/2^(k+1)*sum2;
S
(1)=S
(2);
S
(2)=S(3);
S(3)=T(4)+1/(4-1)*(T(4)-T(3));
C
(1)=C
(2);
C
(2)=S(3)+1/(4^2-1)*(S(3)-S
(2));
R
(2)=C
(2)+1/(4^3-1)*(C
(2)-C
(1));
ifabs(R
(2)-R
(1)) break; end R (1)=R (2); end vpa(R (2),9) &5: 总结matlab自己带有的插值、积分、解线性方程组、解非线性方程及方程组的命令,用实例来展示各命令的用法: &1: Matlab基础应用的实例: 就给出美国1930年到2000年的人口数,根据已有的数据,寻找一种办法来预测一下2000年的美国人口数: 1930 1940 1950 1960 1970 1980 1990 1.23203 1.31669 1.50679 1.79323 2.03212 2.26505 2.49633 解: clc x=1940: 10: 1990; y=[1.232031.316691.506791.793232.032122.265052.49633]; x0=2000; y0=lagrange(x,y,x0) 结果为: y0=3.0394 &2Lagrange插值多项式的程序实例: x=1: 4;y=x.^2; lagrange(x,y,1) ans= 1 lagrange(x,y,2) ans= 4 lagrange(x,y,3) ans= 9 &3Newton插值多项式的程序实例: 例: 给出节点数据 , , , ,作三阶牛顿插值多项式,计算 ,并估计其误差. 解首先将名为newdscg.m的程序保存为M文件,然后在MATLAB工作窗口输入程序 >>symsM,X=[-4,0,1,2];Y=[27,1,2,17];x=-2.345; [y,R,A,C,P]=newdscg(X,Y,x,M) 运行后输出插值y 及其误差限公式R,三阶牛顿插值多项式P及其系数向量C,差商的矩阵A如下 y= 22.3211 R= 1323077530165133/562949953421312*M(即R=2.3503*M) A= 27.0000000 1.0000-6.500000 2.00001.00001.50000 17.000015.00007.00000.9167 C= 0.91674.2500-4.16671.0000 P= 11/12*x^3+17/4*x^2-25/6*x+1 &4龙贝格积分法的上机实现实例: clc clear syms x; n=8;a=0;b=1;R=0.5*10^(-6); T=zeros(n,n); y=sin(x)/x; p=subs(y,x,a); q=subs(y,x,b); if a==0 p=1; end T(1,1)=(b-a)/2*(p+q); for i=2: n f=0; for j=1: 2^(i-2) t=a+((2*j-1)/2^(i-1))*(b-a); z=subs(y,x,t); f=f+z; end T(i,1)=0.5*T(i-1,1)+f*(b-a)/2^(i-1); end for j=2: n n=n-1; for i=1: n T(i,j)=(4^(j-1)*T(i+1,j-1)-T(i,j-1))/(4^(j-1)-1); end end T=vpa(T,7) T=eval(T); n=8; w=ones(1,7); for j=2: n if abs(T(1,j)-T(1,j-1))<=R end end a=T(1,j); a=vpa(a,7) &6 矩阵的直接LU分解法解方程组和微分方程的数值解法上机实现(可以用matlab自带的结法解微分方程组): 编辑一个LU分解函数如下: function[L,U]=Lu(A)%求解线性方程组的三角分解法 %A为方程组的系数矩阵 %L和U为分解后的下三角和上三角矩阵 [n,m]=size(A); ifn~=m error('TherowsandcolumnsofmatrixAmustbeequal! '); return; end%判断矩阵能否LU分解 forii=1: n fori=1: ii forj=1: ii AA(i,j)=A(i,j); end end if(det(AA)==0) error('ThematrixcannotbedividedbyLU! ') return; end end %开始计算,先赋初值 L=eye(n); U=zeros(n,n); %计算U的第一行,L的第一列 fori=1: n U(1,i)=A(1,i); L(i,1)=A(i,1)/U(1,1); end%计算U的第r行,L的第r列 fori=2: n forj=i: n fork=1: i-1 M(k)=L(i,k)*U(k,j); end U(i,j)=A(i,j)-sum(M); end forj=i+1: n fork=1: i-1 M(k)=L(j,k)*U(k,i); end L(j,i)=(A(j,i)-sum(M))/U(i,i); end end 然后,编辑一个通过LU分解法解线性方程组的函数如下 function[L,U,x]=Lu_x(A,d) %三角分解法求解线性方程组,LU法解线性方程组Ax=LUx=d %A为方程组的系数矩阵%d为方程组的右端项 %L和U为分解后的下三角和上三角矩阵 %x为线性方程组的解 [n,m]=size(A); ifn~=m error('TherowsandcolumnsofmatrixAmustbeequal! '); return; end %判断矩阵能否LU分解 forii=1: n fori=1: ii forj=1: ii AA(i,j)=A(i,j); end end if(det(AA)==0) error('ThematrixcannotbedividedbyLU! ') return; end end [L,U]=Lu(A); %直接调用自定义函数,首先将矩阵分解,A=LU %设Ly=d由于L是下三角矩阵,所以可求 y(i)y (1)=d (1); fori=2: n forj=1: i-1 d(i)=d(i)-L(i,j)*y(j); end y(i)=d(i); end %设Ux=y,由于U是上三角矩阵,所以可求x(i)x(n)=y(n)/U(n,n); fori=(n-1): -1: 1 forj=n: -1: i+1 y(i)=y(i)-U(i,j)*x(j); end x(i)=y(i)/U(i,i); end 然后,n=5时,调用自定义函数 [L,U,x]=Lu_x(A,a) 解出: x=0.9999999999896551.0000000000326090.9999999999616251.0000000000199840.999999999996114 [L,U,x]=Lu_x(B,b) 解出: x=0.9999999999999761.0000000000003420.9999999999988191.0000000000014770.999999999999388 微分方程的数值解法上机实现程序: clc u=zeros(11,7); u1=zeros(11,7); pi=3.1415926; h=0.1; k=0.1; r=k/h; fori=1: 11 u(i,1)=(1/8)*sin(pi*(i-1)*h) end fori=1: 11u(i,2)=(1/8)*sin(pi*(i-1)*h)+r*r*(1/8)*sin(pi*(i-1)*h)*(cos(pi*h)-1); end forj=1: 7 u(1,j)=0; end forj=1: 7 u(10,j)=0; end forj=2: 6
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 计算方法 实验 报告