数值分析实习报告

数值分析实习报告

题目：MATLAB算法实习

学生姓名：

指导教师：***

学院：

学号：

2016年12月

目录

一、拉格朗日插值法 (1)

二、最小二乘拟合法 (3)

三、复化辛普森法计算积分 (5)

四、LU分解法 (6)

五、Jacobi迭代法 (8)

六、牛顿迭代法 (9)

七、幂法 (11)

八、欧拉法 (13)

一、拉格朗日插值法

x 0 1 4 9 16

y 0 1 2 3 4

1、建立M文件

function f=language(x,y,x0)

syms t l;

if (length(x)==length(y))

n=length(x);

else

disp('x和y的维数不相等！')

return;

end

h=sym(0);

for i=1:n

l=sym(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;

h=h+l;

end

simplify(h);

if(nargin==3)

f=subs(h,t,x0);

else

f=collect(h);

f=vpa(f,6);

end

2、在命令窗口输入指令并运行:

二、最小二乘拟合法

例题：编制以函数 为基的多项式最小二乘拟合程序，并用于对下列数据作三次多项式最小二乘拟合(取权函数wi≡1)

x 1.0 2.0 3.0 4.0 5.0 6.0 7.0 y

-4.447

-0.452

0.551

0.048

-0.447

0.549

4.552

1、建立主程序文件

x=1.0:1.0:7.0;

y=[-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552]; plot(x,y,'*') xlabel 'x 轴' ylabel 'y 轴' title '散点图' hold on m=6;n=3; A=zeros(n+1); for j=1:n+1 for i=1:n+1 for k=1:m+1

A(j,i)=A(j,i)+x(k)^(j+i-2) end end end;

B=[0 0 0 0]; for j=1:n+1 for i=1:m+1

{}

n

k k

x 0=

B(j)=B(j)+y(i)*x(i)^(j-1)

end

end

B=B';

a=inv(A)*B;

x=[1.0:0.0001:7.0];

z=a(1)+a(2)*x+a(3)*x.^2+a(4)*x.^3;

plot(x,z)

legend('离散点','y=a(1)+a(2)*x+a(3)*x.^2+a(4)*x.^3') title '拟合图'

2、在命令窗口输入指令并运行:

三、复化辛普森法计算积分

例题：利用复化辛普森公式计算积分：xdx

x ln 1

1、建立M 文件

function y=f(x)

y=sqrt(x)*log(x);

2、建立M 文件

function Tn=simp(a,b,n)

h=(b-a)/n;

for k=0:n

x(k+1)=a+k*h; if x(k+1)==0 x(k+1)=10^(-10); end end

T1=h/2*(f(x(1))+f(x(n+1))); for i=2:n

F(i)=h*f(x(i)); end

T2=sum(F); Tn=T1+T2;

3、在命令窗口输入指令并运行:

四、LU 分解法

例题：用LU 分解法解方程组,Ax b =其中

⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=311-2421-17A ，⎥⎥

⎥⎦⎤

⎢⎢⎢⎣⎡=213b

1、建立M 文件

function LU(A,b) %A 为系数矩阵，b 为右端项矩阵 [m,n]=size(A); %初始化矩阵A ，b ，L 和U n=length(b); L=eye(n,n); U=zeros(n,n);

U(1,1:n)=A(1,1:n); %开始进行LU 分解 L(2:n,1)=A(2:n,1)/U(1,1); for k=2:n

U(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);

L(k+1:n,k)=(A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k))/U(k,k); end

L %输出L 矩阵 U %输出U 矩阵

y=zeros(n,1); %开始解方程组Ux=y y(1)=b(1); for k=2:n

y(k)=b(k)-L(k,1:k-1)*y(1:k-1); end x=zeros(n,1); x(n)=y(n)/U(n,n); for k=n-1:-1:1

x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k);