matlab用于计算方法的源程序

1、Newdon迭代法求解非线性方程

function [x k t]=NewdonToEquation(f,df,x0,eps)

%牛顿迭代法解线性方程

%[x k t]=NewdonToEquation(f,df,x0,eps)

%x:近似解

%k:迭代次数

%t:运算时间

%f:原函数，定义为内联函数

?:函数的倒数，定义为内联函数

%x0:初始值

%eps:误差限

%

%应用举例：

%f=inline('x^3+4*x^2-10');

?=inline('3*x^2+8*x');

%x=NewdonToEquation(f,df,1,0.5e-6)

%[x k]=NewdonToEquation(f,df,1,0.5e-6)

%[x k t]=NewdonToEquation(f,df,1,0.5e-6)

%函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquation(f,df,1)

if nargin==3

eps="0".5e-6;

end

tic;

k=0;

while 1

x="x0-f"(x0)./df(x0);

k="k"+1;

if abs(x-x0) < eps || k >30

break;

end

x0=x;

end

t=toc;

if k >= 30

disp('迭代次数太多。');

x="0";

t="0";

end

2、Newdon迭代法求解非线性方程组

function y="NewdonF"(x)

%牛顿迭代法解非线性方程组的测试函数

%定义是必须定义为列向量

y(1,1)=x(1).^2-10*x(1)+x(2).^2+8;

y(2,1)=x(1).*x(2).^2+x(1)-10*x(2)+8;

return;

function y="NewdonDF"(x)

%牛顿迭代法解非线性方程组的测试函数的导数

y(1,1)=2*x(1)-10;

y(1,2)=2*x(2);

y(2,1)=x(2).^+1;

y(2,2)=2*x(1).*x(2)-10;

return;

以上两个函数仅供下面程序的测试

function [x k t]=NewdonToEquations(f,df,x0,eps)

%牛顿迭代法解非线性方程组

%[x k t]=NewdonToEquations(f,df,x0,eps)

%x:近似解

%k:迭代次数

%t:运算时间

%f:方程组（事先定义）

?:方程组的导数（事先定义）

%x0:初始值

%eps:误差限

%

%说明：由于虚参f和df的类型都是函数，使用前需要事先在当前目录下采用函数M文件定义% 另外在使用此函数求解非线性方程组时，需要在函数名前加符号“@”，如下所示

%

%应用举例：

%x0=[0,0];eps=0.5e-6;

%x=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)

%[x k]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)

%[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)

%函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6

%[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)

if nargin==3

eps="0".5e-6;

end

tic;

k=0;

while 1

x="x0-inv"(df(x0))*f(x0);%此处可采用其他方法避免求逆

k="k"+1;

if norm(x-x0) < eps || k > 15

break;

end

x0=x;

end

t=toc;

if k >= 15

disp('迭代次数太多。');

x="zeros"(size(x0));

t="0";

end

3、Lagrange插值法

提供两个程序，采用了不同的方法

function f="InterpLagrange"(x,y,x0)

%构造Lagrange插值多项式

%此函数中借助向量卷积来求Lagrange基函数，运算速度较快

%f=InterpLagrange(x,y,x0)

%f:插值多项式或者是插值多项式在x0处的值

%x:节点

%y:函数值

%x0:某一测试点

%

%调用格式:

%f=InterpLagrange(x,y) 返回插值多项式

%f=InterpLagrange(x,y,x0) 返回插值多项式在点x0处的值

%举例:

%x=[0.32 0.34 0.36];y=[0.314567 0.333487 0.352274];x0=0.33;

%f=InterpLagrange(x,y)

%f=InterpLagrange(x,y,x0)

if length(x)==length(y)

n="length"(x);

else

disp('节点个数和函数值个数不同！')

f=' ';

return;

end

p=0;

for i="1:n"

l="y"(i);

for j="1:n"

if j==i

continue;

end

%利用卷积计算Lagrange基函数

l=conv(l,[1 -x(j)]./(x(i)-x(j)));

end

%p是一向量，表示插值多项式的系数

p="p"+l;

end

if nargin==3

f="polyval"(p,x0);%计算插值多项式在x0处的值

else

f="poly2str"(p,'x');%把插值多项式的向量形式转化为插值多项式的符号形式end

function f="InterpLagrange2"(x,y,x0)

%构造Lagrange插值多项式

%此函数中借助符号运算来求Lagrange基函数，运算速度较慢，不推荐此种方法%f=InterpLagrange2(x,y,x0)

%f:插值多项式或者是插值多项式在x0处的值

%x:节点

%y:函数值

%x0:某一测试点

%

%调用格式:

%f=InterpLagrange2(x,y) 返回插值多项式

%f=InterpLagrange2(x,y,x0) 返回插值多项式在点x0处的值

%举例:

%x=[0.32 0.34 0.36];y=[0.314567 0.333487 0.352274];x0=0.33; %f=InterpLagrange2(x,y)

%f=InterpLagrange2(x,y,x0)

if length(x)==length(y)

n="length"(x);

else

disp('节点个数和函数值个数不同！')

f=' ';

return;

end

syms t;

f=0;

for i="1:n"

l="y"(i);

for j="1:n"

if j==i

continue;

end

l=l*(t-x(j))/(x(i)-x(j));%借助符号运算，计算Lagrange基函数

end

f="f"+l;

simplify(f);%化简多项式

if i==n

if nargin==3

f="subs"(f,'t',x0);%计算插值多项式f在点x0处的值

else

f="collect"(f);%计算插值多项式，展开并合并同类项

f="vpa"(f,6);%设置多项式系数的有效数字

end

end

end

4、Newdon插值法

function f="InterpNewdon"(x,y,x0)

%Newdon插值多项式

%f=InterpNewdon(x,y,x0)

%f:插值多项式或者是插值多项式在x0处的值

%x:节点

%y:函数值

%x0:某一测试点

%

%调用格式

%f=InterpNewdon(x,y) 返回插值多项式

%f=InterpNewdon(x,y,x0) 返回插值多项式在x0点的值

%应用举例:

%x=[1 2 3 4 5];y=[1 4 7 8 6];x0=6;

%f=InterpNewdon(x,y)

%f=InterpNewdon(x,y,x0)

if length(x)==length(y)

n="length"(x);

else

disp('节点个数和函数值个数不同！')

f=' ';

return;

end

A=zeros(n);%初始化差商矩阵

for i="1:n"

A(i,1)=y(i);%差商矩阵的第一列是函数值

end

%计算差商矩阵

%差商矩阵中对角线上的元素为Newdon插值多项式的系数

for j="2:n"

for i="j:n"

A(i,j)=(A(i,j-1)-A(i-1,j-1))/(x(i)-x(i-j+1));

end

end

%求Newdon插值多项式

p=zeros(1,n);

for i="1:n"

p1=A(i,i);%差商矩阵对角线上的元素就是Newdon插值多项式的系数 for j="1:i-1"

p1=conv(p1,[1 -x(j)]);%计算Newdon插值多项式的基项

end

p1=[zeros(1,n-i),p1];%向量相加，维数必须相同。把向量的元素补齐

p="p"+p1;

end

if nargin==3

f="polyval"(p,x0);%计算插值多项式在x0处的值

else

f="poly2str"(p,'x');%把插值多项式的向量形式转化为插值多项式的符号形式end

5、基本Guass消去法求解线性方程组

function x="EqtsBasicGuass"(A,b)

%基本Guass消去法求解线性方程组Ax=b

%x=EqtsBasicGuass(A,b)

%x:解向量，列向量

%A:线性方程组的矩阵

%b:列向量

%

%应用举例:

%A=[2 2 3;4 7 7;-2 4 5]; b=[3;1;-7];

%x=EqtsBasicGuass(A,b)

%检查输入参数

if size(A,1) ~= size(b,1)

disp('输入参数有误！');

x=' ';

return;

end

%(A|b)

A=[A b];

%消去过程

n=size(A,1);

l=zeros(n);

for k="1:n-1"

for i="k"+1:n

l(i,k)=A(i,k)/A(k,k);

end

for i="k"+1:n

for j="k"+1:n+1

A(i,j)=A(i,j)-l(i,k)*A(k,j);

end

for j="1:k"

A(i,j)=0;

end

end

end

%回代过程

x=zeros(n,1);

x(n)=A(n,n+1)/A(n,n);

for i="n-1:-1:1"

y="0";

for j="i"+1:n

y=y+A(i,j)*x(j);

end

x(i)=(A(i,n+1)-y)/A(i,i);

end

return;

6、三角分解法求解线性方程组

function x="EqtsDoolittleLU"(A,b)

%Doolittle分解法求解线性方程组Ax=b

%x=EqtsDoolittleLU(A,b)

%x:解向量，列向量

%A:线性方程组的矩阵

%b:列向量

%

%应用举例:

%A=[6 2 1 -1;2 4 1 0;1 1 4 -1;-1 0 -1 3];b=[6;-1;5;-5]; %x=EqtsDoolittleLU(A,b)

%检查输入参数

if size(A,1) ~= size(b,1)

disp('输入参数有误！');

x=' ';

return;

end

%分解

%把L和U的元素存储在A的相应位置上n=length(b);

for k="1:n"

for j="k:n"

z=0;

for r="1:k-1"

z="z"+A(k,r)*A(r,j);

end

A(k,j)=A(k,j)-z;

end

for i="k"+1:n

z=0;

for r="1:k-1"

z="z"+A(i,r)*A(r,k);

end

A(i,k)=(A(i,k)-z)/A(k,k);

end

end

%求解

x=zeros(size(b));

for i="1:n"

z="0";

for k="1:i-1"

z=z+A(i,k)*x(k);

end

x(i)=b(i)-z;

end

for i="n:-1:1"

z="0";

for k="i"+1:n

z=z+A(i,k)*x(k);

end

x(i)=(x(i)-z)/A(i,i);

end

return

7、追赶法求解三对角线性方程组

function x="EqtsForwardAndBackward"(L,D,U,b)

%追赶法求解三对角线性方程组Ax=b

%x=EqtsForwardAndBackward(L,D,U,b)

%x:三对角线性方程组的解

%L:三对角矩阵的下对角线，行向量

%D:三对角矩阵的对角线，行向量

%U:三对角矩阵的上对角线，行向量

%b:线性方程组Ax=b中的b，列向量

%

%应用举例:

%L=[-1 -2 -3];D=[2 3 4 5];U=[-1 -2 -3];b=[6 1 -2 1]'; %x=EqtsForwardAndBackward(L,D,U,b)

%检查参数的输入是否正确

n=length(D);m=length(b);

n1=length(L);n2=length(U);

if n-n1 ~= 1 || n-n2 ~= 1 || n ~= m

disp('输入参数有误！')

x=' ';

return;

end

%追的过程

for i="2:n"

L(i-1)=L(i-1)/D(i-1);

D(i)=D(i)-L(i-1)*U(i-1);

end

x=zeros(n,1);

x(1)=b(1);

for i="2:n"

x(i)=b(i)-L(i-1)*x(i-1);

end

%赶的过程

x(n)=x(n)/D(n);

for i="n-1:-1:1"

x(i)=(x(i)-U(i)*x(i+1))/D(i);

end

return;

8、主元素的Guass消去法求解线性方程组

function x="EqtsClmnPrimElemGuass"(A,b)

%主元素的Guass消去法求解线性方程组Ax=b

%x=EqtsClmnPrimElemGuass(A,b)

%x:解向量，列向量

%A:线性方程组的矩阵

%b:列向量

%

%应用举例:

%A=[-0.002 2 2;1 0.78125 0;3.996 5.5625 4];

%b=[0.4;1.3816;7.4178];

%x=EqtsClmnPrimElemGuass(A,b)

%检查输入参数

if size(A,1) ~= size(b,1)

disp('输入参数有误！');

x=' ';

return;

end

%(A|b)

A=[A b];

%消去过程

n=size(A,1);

l=zeros(n);

for k="1:n-1"

%换行

[a idx1]=max(abs(A(k:n,k)));%寻找绝对值最大的元素的下标

[b idx2]=min(abs(A(k:n,k)));%寻找绝对值最小的元素的下标 idx1=idx1+k-1;

idx2=idx2+k-1;

for j="1:n"+1

c=A(idx1,j);

A(idx1,j)=A(idx2,j);

A(idx2,j)=c;

end

for i="k"+1:n

l(i,k)=A(i,k)/A(k,k);

end

for i="k"+1:n

for j="k"+1:n+1

A(i,j)=A(i,j)-l(i,k)*A(k,j);

end

for j="1:k"

A(i,j)=0;

end

end

end

%回代过程

x=zeros(n,1);

x(n)=A(n,n+1)/A(n,n);

for i="n-1:-1:1"

y="0";

for j="i"+1:n

y=y+A(i,j)*x(j);

end

x(i)=(A(i,n+1)-y)/A(i,i);

end

return;

9、Euler法求解常微分方程

function outXY="ODEEuler"(f,x0,y0,h,PointNum)

%简单欧拉法求解常微分方程dy/dx=f

%outXY=ODEEuler(f,x0,y0,h,PointNum)

%outXY:所取点的横纵坐标。第一列为横坐标，第二列为纵坐标

%f:函数f(x,y)，可利用脚本函数文件事先定义，也可利用内联函数%x0:初始值的横坐标

%y0:初始值的纵坐标

%h:步长

%PointNum:计算步数，默认为30

%

%应用举例:

%f=inline('y-2*x/y','x','y');

%out=ODEEuler(f,0,1,0.1,10)

if nargin==4

PointNum="30";

end

outXY=zeros(PointNum+1,2);%初始化

outXY(1,1)=x0;

outXY(1,2)=y0;

for i="1:PointNum"

outXY(i+1,2)=outXY(i,2)+h*f(outXY(i,1),outXY(i,2));%简单Euler公式 outXY(i+1,1)=outXY(i,1)+h;

end

10、二分法解非线性方程

function [x k t]=DichotomyToEquation(f,a,b,eps)

%使用二分法解非线性方程

%[x k t]=DichotomyToEquation(f,a,b,eps)

%x:近似解

%k:二分次数

%t:运算时间

%f:函数，定义为内联函数

%a,b:区间端点

%eps:误差限

%

%应用举例:

%f=inline('x^3+4*x^2-10');

%x=DichotomyToEquation(f,1,2,0.5e-6)

%[x k]=DichotomyToEquation(f,1,2,0.5e-6)

%[x k t]=DichotomyToEquation(f,1,2,0.5e-6)

%函数的最后一个参数也可以不写，默认情况下，eps=0.5e-6

%[x k t]=DichotomyToEquation(f,1,2)

if nargin==3

eps="0".5e-6;

end

tic;

if f(a)*f(b) > 0

disp('区间太大或在此区间内无零点。');

else

k="0";%记录二分次数

x=(b+a)/2;

k=k+1;

if abs(x-a) < eps || k > 30

break;

end

if f(a)*f(x) < 0

b="x";

else

a="x";

end

end

t="toc";

x=(b+a)/2;

if k >= 30

disp('迭代次数太多。');

x=0;

t=0;

end

end

11、弦割法求解非线性方程

function [x k t]=ChordsecantToEquation(f,x0,x1,eps) %弦割法求解非线性方程

%[x k t]=ChordsecantToEquation(f,x0,x1,eps)

%x:近似解

%k:迭代次数

%t:运算时间

%f:原函数，定义为内联函数

%x0,x1:初始值

%eps:误差限

%

%应用举例：

%f=inline('x^3+4*x^2-10');

%x=ChordsecantToEquation(f,1,2,0.5e-6)

%[x k]=ChordsecantToEquation(f,1,2,0.5e-6)

%[x k t]=ChordsecantToEquation(f,1,2,0.5e-6)

%函数的最后一个参数也可以不写，默认情况下，eps=0.5e-6 %[x k t]=ChordsecantToEquation(f,1,2)

eps="0".5e-6;

end

tic;

k=0;

while 1

x="x1-f"(x1)*(x1-x0)./(f(x1)-f(x0));

k="k"+1;

if abs(x-x1) < eps || k > 30

break;

end

x0=x1;

x1=x;

end

t=toc;

if k >= 30

disp('迭代次数太多。');

x="0";

t="0";

end

function y="NewdonDampingDF"(x)

%带阻尼因子的牛顿迭代法测试，方程组的Jacobi矩阵y(1,1)=2*x(1)-10;

y(1,2)=2*x(2);

y(2,1)=x(2)^2+1;

y(2,2)=2*x(1)*x(2)-10;

return;

以上两个函数用于下列函数的测试

function [x k t]=NewdonDampingToEquations(f,df,x0,yita,eps)

%带阻尼因子的Newdon迭代格式求解非线性方程组

%[x k t]=NewdonDampingToEquations(f,df,x0,yita,eps)

%x:近似解

%k:迭代次数

%t:运算时间

%f:方程组

?:方程组的Jacobi矩阵

%x0:初始值

%yita:阻尼因子，为了克服Jacobi矩阵的奇异性

%eps:误差限

%

%应用举例：

%x0=[2.5;2.5];yita=1e-5;eps=0.5e-6;

%x=NewdonDampingToEquations(@NewdonDampingF,@NewdonDampingDF,x0,yita,eps) %[x

k]=NewdonDampingToEquations(@NewdonDampingF,@NewdonDampingDF,x0,yita,eps) %[x k

t]=NewdonDampingToEquations(@NewdonDampingF,@NewdonDampingDF,x0,yita,eps) %函数的最后两个参数也可以不写，默认情况下，yita=1e-4;eps=0.5e-6

%[x k t]=NewdonDampingToEquations(@NewdonDampingF,@NewdonDampingDF,x0)

%[x k t]=NewdonDampingToEquations(@NewdonDampingF,@NewdonDampingDF,x0,yita)

if nargin==3

yita="1e-4";

eps="0".5e-6;

end

if nargin==4

eps="0".5e-6;

end

I=eye(size(df(x0)));

tic;

k=0;

while 1

x="x0-inv"(df(x0)+yita*I)*f(x0);%此处可采用其他方法避免求逆

k="k"+1;

if norm(x-x0) < eps || k > 30

break;

end

x0=x;

end

t=toc;

if k >= 30

disp('迭代次数太多。');

x="zeros"(size(x0));

t="0";

end

13、牛顿下山法求解非线性方程组

测试函数见2（Newdon迭代法求解非线性方程组）

function [x k t]=NewdonDescendToEquations(f,df,x0,omiga,eps)

%下降牛顿迭代法求解非线性方程组

%[x k t]=NewdonDescendToEquations(f,df,x0,eps)

%x:近似解

%k:迭代次数

%t:运算时间

%f:方程组

?:方程组的Jacobi矩阵

%x0:初始值

%omiga:下降因子，通常在区间(0,1)内选择。当omiga=1时，即为Newdon迭代格式%eps:误差限

%

%应用举例：

%x0=[0;0];eps=0.5e-6;

%x=NewdonDescendToEquations(@NewdonF,@NewdonDF,x0,omiga,eps)

%[x k]=NewdonDescendToEquations(@NewdonF,@NewdonDF,x0,omiga,eps) %[x k t]=NewdonDescendToEquations(@NewdonF,@NewdonDF,x0,omiga,eps) %函数的最后两个参数也可以不写，默认情况下，omiga=0.8;eps=0.5e-6

%[x k t]=NewdonDescendToEquations(@NewdonF,@NewdonDF,x0)

%[x k t]=NewdonDescendToEquations(@NewdonF,@NewdonDF,x0,omiga)

if nargin==3

omiga="0".8;

eps="0".5e-6;

end

if nargin==4

eps="0".5e-6;

end

tic;

k=0;

while 1

x="x0-omiga"*inv(df(x0))*f(x0);%此处可采用其他方法避免求逆

k="k"+1;

if norm(x-x0) < eps || k > 30

break;

end

x0=x;

end

t=toc;

if k >= 30

disp('迭代次数太多。');

x="zeros"(size(x0));

t="0";

end

14、简化牛顿法求解非线性方程组

测试函数见2（Newdon迭代法求解非线性方程组）

function [x k t]=NewdonSimplifyToEquations(f,df,x0,eps)

%简化牛顿格式求解非线性方程组

%[x k t]=NewdonSimplifyToEquations(f,df,x0,eps)

%x:近似解

%k:迭代次数

%t:运算时间

%f:方程组

?:方程组的Jacobi矩阵

%x0:初始值

%eps:误差限

%

%应用举例：

%x0=[0;0];eps=0.5e-6;

%x=NewdonSimplifyToEquations(@NewdonF,@NewdonDF,x0,eps)

%[x k]=NewdonSimplifyToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k t]=NewdonSimplifyToEquations(@NewdonF,@NewdonDF,x0,eps) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6

%[x k t]=NewdonSimplifyToEquations(@NewdonF,@NewdonDF,x0)

if nargin==3

eps="0".5e-6;

end

x_const=x0;

tic;

k=0;

A=inv(df(x_const));

while 1

x="x0-A"*f(x0);%此处可采用其他方法避免求逆

k="k"+1;

if norm(x-x0) < eps || k > 30

break;

end

x0=x;

end

t=toc;

if k >= 30

disp('迭代次数太多。');

x="zeros"(size(x0));

t="0";

end

15、逆Broyden秩1方法求解非线性方程组

测试函数见2（Newdon迭代法求解非线性方程组）

function [x k t]=Broyden1InvToEquations(f,df,x0,eps)

%逆Broyden秩1方法求解非线性方程组

%function [x k t]=Broyden1InvToEquations(f,df,x0,eps)

%x:近似解

%k:迭代次数

%t:运算时间

%f:方程组

?:方程组的Jacobi矩阵

%x0:初始值

%eps:误差限

%

%应用举例：

%x0=[0;0];eps=0.5e-6;

%x=Broyden1InvToEquations(@NewdonF,@NewdonDF,x0,eps)

%[x k]=Broyden1InvToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k t]=Broyden1InvToEquations(@NewdonF,@NewdonDF,x0,eps) %函数的最后一个参数也可以不写，默认情况下，eps=0.5e-6

%[x k t]=Broyden1InvToEquations(@NewdonF,@NewdonDF,x0)

if nargin==3

eps="0".5e-6;

end

tic;

k=0;

B0=inv(df(x0));

while 1

x="x0-B0"*f(x0);

k="k"+1;

if norm(x-x0) < eps || k> 30

break;

end

s="x-x0";

y="f"(x)-f(x0);

B="B0"+(s-B0*y)*s'*B0/(s'*B0*y);

x0=x;

B0=B;

end

t=toc;

if k >= 30

disp('迭代次数太多。');

x="zeros"(size(x0));

t="0";

end

16、Jacobi迭代法求解线性方程组

function [x k]=EqtsJacobi(A,b,x0,eps)

%Jacobi迭代法求解线性方程组Ax=b

%[x k]=EqtsJacobi(A,b,x0,eps)

%x:解向量，列向量

%k:迭代次数

%A:系数矩阵

%b:列向量

matlab源代码实例

1.硬币模拟试验 源代码： clear; clc; head_count=0; p1_hist= [0]; p2_hist= [0]; n = 1000; p1 = 0.3; p2=0.03; head = figure(1); rand('seed',sum(100*clock)); fori = 1:n tmp = rand(1); if(tmp<= p1) head_count = head_count + 1; end p1_hist (i) = head_count /i; end figure(head); subplot(2,1,1); plot(p1_hist); grid on; hold on; xlabel('重复试验次数'); ylabel('正面向上的比率'); title('p=0.3试验次数N与正面向上比率的函数图'); head_count=0; fori = 1:n tmp = rand(1); if(tmp<= p2) head_count = head_count + 1; end p2_hist (i) = head_count /i; end figure(head); subplot(2,1,2); plot(p2_hist); grid on; hold on; xlabel('重复试验次数'); ylabel('正面向上的比率'); title('p=0.03试验次数N与正面向上比率的函数图'); 实验结果：

2.不同次数的随机试验均值方差比较 源代码： clear ; clc; close; rand('seed',sum(100*clock)); Titles = ['n=5时' 'n=20时' 'n=25时' 'n=50时' 'n=100时']; Titlestr = cellstr(Titles); X_n_bar=[0]; %the samples of the X_n_bar X_n=[0]; %the samples of X_n N=[5,10,25,50,100]; j=1; num_X_n = 100; num_X_n_bar = 100; h_X_n_bar = figure(1);

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

6-1 试验目的计算特征值，实现算法 试验容：随机产生一个10阶整数矩阵，各数均在-5和5之间。 （1） 用MATLAB 函数“eig ”求矩阵全部特征值。 （2） 用幂法求A 的主特征值及对应的特征向量。 （3） 用基本QR 算法求全部特征值（可用MATLAB 函数“qr ”实现矩阵的QR 分解）。 原理 幂法：设矩阵A 的特征值为12n ||>||||λλλ≥???≥并设A 有完全的特征向量系12,,,n χχχ???(它们线性无关)，则对任意一个非零向量0n V R ∈所构造的向量序列1k k V AV -=有11()lim ()k j k k j V V λ→∞ -=， 其中()k j V 表示向量的第j 个分量。 为避免逐次迭代向量k V 不为零的分量变得很大（1||1λ>时）或很小（1||1λ<时），将每一步的k V 按其模最大的元素进行归一化。具体过程如下： 选择初始向量0V ，令1max(),,,1k k k k k k k V m V U V AU k m +===≥，当k 充分大时1111,max()max() k k U V χλχ+≈ ≈。 QR 法求全部特征值： 111 11222 111 ,1,2,3,k k k k k A A Q R R Q A Q R k R Q A Q R +++==????==??=???? ??????==?? 由于此题的矩阵是10阶的，上述算法计算时间过长，考虑采用改进算法——移位加速。迭 代格式如下： 1 k k k k k k k k A q I Q R A R Q q I +-=?? =+? 计算k A 右下角的二阶矩阵() () 1,1 1,() (),1 ,k k n n n n k k n n n n a a a a ----?? ? ??? 的特征值()()1,k k n n λλ-，当()()1,k k n n λλ-为实数时，选k q 为()()1,k k n n λλ-中最接近(),k n n a 的。 程序

matlab语音识别系统(源代码)最新版

matlab语音识别系统(源代码)最新版

目录 一、设计任务及要求 (1) 二、语音识别的简单介绍 2.1语者识别的概念 (2) 2.2特征参数的提取 (3) 2.3用矢量量化聚类法生成码本 (3) 2.4VQ的说话人识别 (4) 三、算法程序分析 3.1函数关系 (4) 3.2代码说明 (5) 3.2.1函数mfcc (5) 3.2.2函数disteu (5) 3.2.3函数vqlbg (6) 3.2.4函数test (6) 3.2.5函数testDB (7) 3.2.6 函数train (8) 3.2.7函数melfb (8) 四、演示分析 (9) 五、心得体会 (11) 附：GUI程序代码 (12)

一、设计任务及要求 用MATLAB实现简单的语音识别功能； 具体设计要求如下： 用MATLAB实现简单的数字1~9的语音识别功能。 二、语音识别的简单介绍 基于VQ的说话人识别系统，矢量量化起着双重作用。在训练阶段，把每一个说话者所提取的特征参数进行分类，产生不同码字所组成的码本。在识别(匹配)阶段，我们用VQ方法计算平均失真测度(本系统在计算距离d时，采用欧氏距离测度)，从而判断说话人是谁。 语音识别系统结构框图如图1所示。 图1 语音识别系统结构框图 2.1语者识别的概念 语者识别就是根据说话人的语音信号来判别说话人的身份。语音是人的自然属性之一，由于说话人发音器官的生理差异以及后天形成的行为差异，每个人的语音都带有强烈的个人色彩，这就使得通过分析语音信号来识别说话人成为可能。用语音来鉴别说话人的身份有着许多独特的优点，如语音是人的固有的特征，不会丢失或遗忘；语音信号的采集方便，系统设备成本低；利用电话网络还可实现远程客户服务等。因此，近几年来，说话人识别越来越多的受到人们的重视。与其他生物识别技术如指纹识别、手形识别等相比较，说话人识别不仅使用方便，而且属于非接触性，容易被用户接受，并且在已有的各种生物特征识别技术中，是唯一可以用作远程验证的识别技术。因此，说话人识别的应用前景非常广泛：今天，说话人识别技术已经关系到多学科的研究领域，不同领域中的进步都对说话人识别的发展做出了贡献。说话人识别技术是集声学、语言学、计算机、信息处理和人工智能等诸多领域的一项综合技术，应用需求将十分广阔。在吃力语音信号的时候如何提取信号中关键的成分尤为重要。语音信号的特征参数的好坏直接导致了辨别的准确性。

计算方法_全主元消去法_matlab程序

%求四阶线性方程组的MA TLAB程序 clear Ab=[0.001 2 1 5 1; 3 - 4 0.1 -2 2; 2 -1 2 0.01 3; 1.1 6 2.3 9 4];%增广矩阵 num=[1 2 3 4];%未知量x的对应序号 for i=1:3 A=abs(Ab(i:4,i:4));%系数矩阵取绝对值 [r,c]=find(A==max(A(:))); r=r+i-1;%最大值对应行号 c=c+i-1;%最大值对应列号 q=Ab(r,:),Ab(r,:)=Ab(i,:),Ab(i,:)=q;%行变换 w=Ab(:,c),Ab(:,c)=Ab(:,i),Ab(:,i)=w;%列变换 n=num(i),num(i)=num(c),num(c)=n;%列变换引起未知量x次序变化for j=i:3 Ab(j+1,:)=-Ab(j+1,i)*Ab(i,:)/Ab(i,i)+Ab(j+1,:);%消去过程 end end %最后得到系数矩阵为上三角矩阵 %回代算法求解上三角形方程组 x(4)=Ab(4,5)/Ab(4,4); x(3)=(Ab(3,5)-Ab(3,4)*x(4))/Ab(3,3); x(2)=(Ab(2,5)-Ab(2,3)*x(3)-Ab(2,4)*x(4))/Ab(2,2); x(1)=(Ab(1,5)-Ab(1,2)*x(2)-Ab(1,3)*x(3)-Ab(1,4)*x(4))/Ab(1,1); for s=1:4 fprintf('未知量x%g =%g\n',num(s),x(s)) end %验证如下 %A=[0.001 2 1 5 1; 3 -4 0.1 -2 2;2 -1 2 0.01 3; 1.1 6 2.3 9 4]; %b=[1 2 3 4]'; %x=A\b; %x1= 1.0308 %x2= 0.3144 %x3= 0.6267 %x4= -0.0513

基于MATLAB的潮流计算源程序代码(优.选)

%*************************电力系统直角坐标系下的牛顿拉夫逊法潮流计算********** clear clc load E:\data\IEEE014_Node.txt Node=IEEE014_Node; weishu=size(Node); nnum=weishu(1,1); %节点总数 load E:\data\IEEE014_Branch.txt branch=IEEE014_Branch; bwei=size(branch); bnum=bwei(1,1); %支路总数 Y=(zeros(nnum)); Sj=100; %********************************节点导纳矩阵******************************* for m=1:bnum; s=branch(m,1); %首节点 e=branch(m,2); %末节点 R=branch(m,3); %支路电阻 X=branch(m,4); %支路电抗 B=branch(m,5); %支路对地电纳 k=branch(m,6); if k==0 %无变压器支路情形 Y(s,e)=-1/(R+j*X); %互导纳 Y(e,s)=Y(s,e); end if k~=0 %有变压器支路情形 Y(s,e)=-(1/((R+j*X)*k)); Y(e,s)=Y(s,e); Y(s,s)=-(1-k)/((R+j*X)*k^2); Y(e,e)=-(k-1)/((R+j*X)*k); %对地导纳 end Y(s,s)=Y(s,s)-j*B/2; Y(e,e)=Y(e,e)-j*B/2; %自导纳的计算情形 end for t=1:nnum; Y(t,t)=-sum(Y(t,:))+Node(t,12)+j*Node(t,13); %求支路自导纳 end G=real(Y); %电导 B=imag(Y); %电纳 %******************节点分类************************************* * pq=0; pv=0; blancenode=0; pqnode=zeros(1,nnum); pvnode=zeros(1,nnum); for m=1:nnum; if Node(m,2)==3 blancenode=m; %平衡节点编号 else if Node(m,2)==0 pq=pq+1; pqnode(1,pq)=m; %PQ 节点编号 else if Node(m,2)==2 pv=pv+1; pvnode(1,pv)=m; %PV 节点编号 end end end end %*****************************设置电压初值********************************** Uoriginal=zeros(1,nnum); %对各节点电压矩阵初始化 for n=1:nnum Uoriginal(1,n)=Node(n,9); %对各点电压赋初值 if Node(n,9)==0;

最常用的matlab图像处理的源代码

最常用的一些图像处理Matlab源代 码 #1：数字图像矩阵数据的显示及其傅立叶变换 #2：二维离散余弦变换的图像压缩 #3：采用灰度变换的方法增强图像的对比度 #4：直方图均匀化 #5：模拟图像受高斯白噪声和椒盐噪声的影响 #6：采用二维中值滤波函数medfilt2对受椒盐噪声干扰的图像滤波 #7：采用MATLAB中的函数filter2对受噪声干扰的图像进行均值滤波 #8：图像的自适应魏纳滤波 #9：运用5种不同的梯度增强法进行图像锐化 #10：图像的高通滤波和掩模处理 #11：利用巴特沃斯（Butterworth）低通滤波器对受噪声干扰的图像进行平滑处理 #12：利用巴特沃斯（Butterworth）高通滤波器对受噪声干扰的图像进行平滑处理 1.数字图像矩阵数据的显示及其傅立叶变换 f=zeros(30,30); f(5:24,13:17)=1; imshow(f, 'notruesize'); F=fft2(f,256,256); % 快速傅立叶变换算法只能处矩阵维数为2的幂次，f矩阵不 % 是，通过对f矩阵进行零填充来调整 F2=fftshift(F); % 一般在计算图形函数的傅立叶变换时，坐标原点在 % 函数图形的中心位置处，而计算机在对图像执行傅立叶变换 % 时是以图像的左上角为坐标原点。所以使用函数fftshift进 %行修正，使变换后的直流分量位于图形的中心； figure,imshow(log(abs(F2)),[-1 5],'notruesize');

2 二维离散余弦变换的图像压缩I=imread('cameraman.tif'); % MATLAB自带的图像imshow(I); clear;close all I=imread('cameraman.tif'); imshow(I); I=im2double(I); T=dctmtx(8); B=blkproc(I,[8 8], 'P1*x*P2',T,T'); Mask=[1 1 1 1 0 0 0 0 1 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]; B2=blkproc(B,[8 8],'P1.*x',Mask); % 此处为点乘(.*) I2=blkproc(B2,[8 8], 'P1*x*P2',T',T); figure,imshow(I2); % 重建后的图像 3.采用灰度变换的方法增强图像的对比度I=imread('rice.tif'); imshow(I); figure,imhist(I); J=imadjust(I,[0.15 0.9], [0 1]); figure,imshow(J); figure,imhist(J);

王能超 计算方法——算法设计及MATLAB实现课后代码

第一章插值方法 1.1Lagrange插值 1.2逐步插值 1.3分段三次Hermite插值 1.4分段三次样条插值 第二章数值积分 2.1 Simpson公式 2.2 变步长梯形法 2.3 Romberg加速算法 2.4 三点Gauss公式 第三章常微分方程德差分方法 3.1 改进的Euler方法 3.2 四阶Runge-Kutta方法 3.3 二阶Adams预报校正系统 3.4 改进的四阶Adams预报校正系统 第四章方程求根 4.1 二分法 4.2 开方法 4.3 Newton下山法 4.4 快速弦截法 第五章线性方程组的迭代法 5.1 Jacobi迭代 5.2 Gauss-Seidel迭代 5.3 超松弛迭代 5.4 对称超松弛迭代 第六章线性方程组的直接法 6.1 追赶法 6.2 Cholesky方法 6.3 矩阵分解方法 6.4 Gauss列主元消去法

第一章插值方法 1.1Lagrange插值 计算Lagrange插值多项式在x=x0处的值. MATLAB文件：（文件名：Lagrange_eval.m）function [y0,N]= Lagrange_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Lagrange插值多项式在x0处的值 %N是Lagrange插值函数的权系数 m=length(X); N=zeros(m,1); y0=0; for i=1:m N(i)=1; for j=1:m if j~=i; N(i)=N(i)*(x0-X(j))/(X(i)-X(j)); end end y0=y0+Y(i)*N(i); end 用法》X=[…];Y=[…]; 》x0= ; 》[y0,N]= Lagrange_eval(X,Y,x0) 1.2逐步插值 计算逐步插值多项式在x=x0处的值. MATLAB文件：（文件名：Neville_eval.m）function y0=Neville_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Neville逐步插值多项式在x0处的值 m=length(X); P=zeros(m,1); P1=zeros(m,1); P=Y; for i=1:m P1=P; k=1; for j=i+1:m k=k+1;

BP神经网络matlab源程序代码

close all clear echo on clc % NEWFF——生成一个新的前向神经网络 % TRAIN——对 BP 神经网络进行训练 % SIM——对 BP 神经网络进行仿真 % 定义训练样本 % P为输入矢量 P=[0.7317 0.6790 0.5710 0.5673 0.5948;0.6790 0.5710 0.5673 0.5948 0.6292; ... 0.5710 0.5673 0.5948 0.6292 0.6488;0.5673 0.5948 0.6292 0.6488 0.6130; ... 0.5948 0.6292 0.6488 0.6130 0.5654; 0.6292 0.6488 0.6130 0.5654 0.5567; ... 0.6488 0.6130 0.5654 0.5567 0.5673;0.6130 0.5654 0.5567 0.5673 0.5976; ... 0.5654 0.5567 0.5673 0.5976 0.6269;0.5567 0.5673 0.5976 0.6269 0.6274; ... 0.5673 0.5976 0.6269 0.6274 0.6301;0.5976 0.6269 0.6274 0.6301 0.5803; ... 0.6269 0.6274 0.6301 0.5803 0.6668;0.6274 0.6301 0.5803 0.6668 0.6896; ... 0.6301 0.5803 0.6668 0.6896 0.7497]; % T为目标矢量 T=[0.6292 0.6488 0.6130 0.5654 0.5567 0.5673 0.5976 ... 0.6269 0.6274 0.6301 0.5803 0.6668 0.6896 0.7497 0.8094]; % Ptest为测试输入矢量 Ptest=[0.5803 0.6668 0.6896 0.7497 0.8094;0.6668 0.6896 0.7497 0.8094 0.8722; ... 0.6896 0.7497 0.8094 0.8722 0.9096]; % Ttest为测试目标矢量 Ttest=[0.8722 0.9096 1.0000]; % 创建一个新的前向神经网络 net=newff(minmax(P'),[12,1],{'logsig','purelin'},'traingdm'); % 设置训练参数 net.trainParam.show = 50; net.trainParam.lr = 0.05; net.trainParam.mc = 0.9; net.trainParam.epochs = 5000; net.trainParam.goal = 0.001; % 调用TRAINGDM算法训练 BP 网络 [net,tr]=train(net,P',T); % 对BP网络进行仿真 A=sim(net,P'); figure; plot((1993:2007),T,'-*',(1993:2007),A,'-o'); title('网络的实际输出和仿真输出结果，*为真实值，o为预测值'); xlabel('年份'); ylabel('客运量'); % 对BP网络进行测试 A1=sim(net,Ptest');

(整理)matlab16常用计算方法.

常用计算方法 1．超越方程的求解 一超越方程为 x (2ln x – 3) -100 = 0 求超越方程的解。 [算法]方法一：用迭代算法。将方程改为 01002ln()3 x x =- 其中x 0是一个初始值，由此计算终值x 。取最大误差为e = 10-4，当| x - x 0| > e 时，就用x 的值换成x 0的值，重新进行计算；否则| x - x 0| < e 为止。 [程序]P1_1abs.m 如下。 %超越方程的迭代算法 clear %清除变量 x0=30; %初始值 xx=[]; %空向量 while 1 %无限循环 x=100/(2*log(x0)-3); %迭代运算 xx=[xx,x]; %连接结果 if length(xx)>1000,break ,end %如果项数太多则退出循环（暗示发散） if abs(x0-x)<1e-4,break ,end %当精度足够高时退出循环 x0=x; %替换初值 end %结束循环 figure %创建图形窗口 plot(xx,'.-','LineWidth',2,'MarkerSize',12)%画迭代线'.-'表示每个点用.来表示，再用线连接 grid on %加网格 fs=16; %字体大小 title('超越方程的迭代折线','fontsize',fs)%标题 xlabel('\itn','fontsize',fs) %x 标签 ylabel('\itx','fontsize',fs) %y 标签 text(length(xx),xx(end),num2str(xx(end)),'fontsize',fs)%显示结果 [图示]用下标作为自变量画迭代的折线。如P0_20_1图所示，当最大误差为10-4时，需要迭代19次才能达到精度，超越方程的解为27.539。 [算法]方法二：用求零函数和求解函数。将方程改为函数 100()2ln()3f x x x =-- MATLAB 求零函数为fzero ，fzero 函数的格式之一是 x = fzero(f,x0) 其中，f 表示求解的函数文件，x0是估计值。fzero 函数的格式之二是 x = fzero(f,[x1,x2])

Matlab源程序代码

正弦波的源程序： (一)，用到的函数 1,f2t函数 function x=f2t(X) global dt df t f T N %x=f2t(X) %x为时域的取样值矢量 %X为x的傅氏变换 %X与x长度相同并为2的整幂 %本函数需要一个全局变量dt(时域取样间隔) X=[X(N/2+1:N),X(1:N/2)]; x=ifft(X)/dt; end 2,t2f函数。 function X=t2f(x) global dt df N t f T %X=t2f(x) %x为时域的取样值矢量 %X为x的傅氏变换 %X与x长度相同,并为2的整幂。 %本函数需要一个全局变量dt(时域取样间隔) H=fft(x); X=[H(N/2+1:N),H(1:N/2)]*dt; end (二),主程序。 1，%（1）绘出正弦信号波形及频谱 global dt df t f N close all k=input('取样点数＝2^k, k取10左右'); if isempty(k), k=10; end f0=input('f0=取1(kz)左右'); if isempty(f0), f0=1; end N=2^k; dt=0.01; %ms df=1/(N*dt); %KHz T=N*dt; %截短时间

Bs=N*df/2; %系统带宽 f=[-Bs+df/2:df:Bs]; %频域横坐标 t=[-T/2+dt/2:dt:T/2]; %时域横坐标 s=sin(2*pi*f0*t); %输入的正弦信号 S=t2f(s); %S是s的傅氏变换 a=f2t(S); %a是S的傅氏反变换 a=real(a); as=abs(S); subplot(2,1,1) %输出的频谱 plot(f,as,'b'); grid axis([-2*f0,+2*f0,min(as),max(as)]) xlabel('f (KHz)') ylabel('|S(f)| (V/KHz)') %figure(2) subplot(2,1,2) plot(t,a,'black') %输出信号波形画图grid axis([-2/f0,+2/f0,-1.5,1.5]) xlabel('t(ms)') ylabel('a(t)(V)') gtext('频谱图') 最佳基带系统的源程序： (一)，用到的函数 f2t函数和t2f函数。代码>> (二),主程序 globaldt t f df N T close all clear Eb_N0 Pe k=input('取样点数＝2^k, k取13左右'); if isempty(k), k=13; end z=input('每个信号取样点数＝2^z, z

matlab用于计算方法的源程序

1、Newdon迭代法求解非线性方程 function [x k t]=NewdonToEquation(f,df,x0,eps) %牛顿迭代法解线性方程 %[x k t]=NewdonToEquation(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:原函数，定义为内联函数 ?:函数的倒数，定义为内联函数 %x0:初始值 %eps:误差限 % %应用举例： %f=inline('x^3+4*x^2-10'); ?=inline('3*x^2+8*x'); %x=NewdonToEquation(f,df,1,0.5e-6) %[x k]=NewdonToEquation(f,df,1,0.5e-6) %[x k t]=NewdonToEquation(f,df,1,0.5e-6) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquation(f,df,1) if nargin==3 eps="0".5e-6; end tic; k=0; while 1 x="x0-f"(x0)./df(x0); k="k"+1; if abs(x-x0) < eps || k >30 break; end x0=x; end t=toc; if k >= 30 disp('迭代次数太多。'); x="0"; t="0"; end

2、Newdon迭代法求解非线性方程组 function y="NewdonF"(x) %牛顿迭代法解非线性方程组的测试函数 %定义是必须定义为列向量 y(1,1)=x(1).^2-10*x(1)+x(2).^2+8; y(2,1)=x(1).*x(2).^2+x(1)-10*x(2)+8; return; function y="NewdonDF"(x) %牛顿迭代法解非线性方程组的测试函数的导数 y(1,1)=2*x(1)-10; y(1,2)=2*x(2); y(2,1)=x(2).^+1; y(2,2)=2*x(1).*x(2)-10; return; 以上两个函数仅供下面程序的测试 function [x k t]=NewdonToEquations(f,df,x0,eps) %牛顿迭代法解非线性方程组 %[x k t]=NewdonToEquations(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:方程组（事先定义） ?:方程组的导数（事先定义） %x0:初始值 %eps:误差限 % %说明：由于虚参f和df的类型都是函数，使用前需要事先在当前目录下采用函数M文件定义% 另外在使用此函数求解非线性方程组时，需要在函数名前加符号“@”，如下所示 % %应用举例： %x0=[0,0];eps=0.5e-6; %x=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)

0计算方法及MATLAB实现简明讲义课件PPS8-1欧拉龙格法

第8章 常微分方程初值问题数值解法 8.1 引言 8.2 欧拉方法 8.3 龙格-库塔方法 8.4 单步法的收敛性与稳定性 8.5 线性多步法

8.1 引 言 考虑一阶常微分方程的初值问题 00(,),[,],(). y f x y x a b y x y '=∈=（1.1） （1.2） 如果存在实数 ，使得 121212(,)(,).,R f x y f x y L y y y y -≤-?∈（1.3） 则称 关于 满足李普希茨(Lipschitz ）条件， 称为 的李普希茨常数（简称Lips.常数）. 0>L f y L f (参阅教材386页)

计算方法及MATLAB 实现 所谓数值解法，就是寻求解 在一系列离散节点 )(x y <<<<<+121n n x x x x 上的近似值 . ,,,,,121+n n y y y y 相邻两个节点的间距 称为步长. n n n x x h -=+1 如不特别说明，总是假定 为定数， ),2,1( ==i h h i 这时节点为 . ) ,2,1,0(0 =+=i nh x x n 初值问题(1.1),(1.2)的数值解法的基本特点是采取 “步进式”. 即求解过程顺着节点排列的次序一步一步地向前推进. 00(,),[,], (). y f x y x a b y x y '=∈=

描述这类算法，只要给出用已知信息 ,,,21--n n n y y y 计算 的递推公式. 1+n y 一类是计算 时只用到前一点的值 ，称为单步法. 1+n y n y 另一类是用到 前面 点的值 ， 1+n y k 11,,,+--k n n n y y y 称为 步法. k 其次，要研究公式的局部截断误差和阶，数值解 与 精确解 的误差估计及收敛性，还有递推公式的计算 稳定性等问题. n y )(n x y 首先对方程 离散化，建立求数值解的递推 公式. ),(y x f y ='

主成分分析matlab源程序代码

263.862 1.61144 2.754680.266575 268.764 2.07218 2.617560.182597 261.196 1.59769 2.350370.182114 248.708 2.09609 2.852790.257724 253.365 1.69457 2.94920.189702 268.434 1.56819 2.781130.13252 258.741 2.14653 2.691110.136469 244.192 2.02156 2.226070.298066 219.738 1.61224 1.885990.166298 244.702 1.91477 2.259450.187569 245.286 2.12499 2.352820.161602 251.96 1.83714 2.535190.240271 251.164 1.74167 2.629610.211887 251.824 2.00133 2.626650.211991 257.68 2.14878 2.656860.203846] stdr=std(dataset);%求个变量的标准差 [n,m]=size(dataset);%定义矩阵行列数 sddata=dataset./stdr(ones(n,1),:);%将原始数据采集标准化 sddata%输出标准化数据 [p,princ,eigenvalue,t2]=princomp(sddata);%调用前三个主成分系数 p3=p(:,1:3);%提取前三个主成分得分系数,通过看行可以看出对应的原始数据的列，每个列在每个主成分的得分 p3%输出前三个主成分得分系数 sc=princ(:,1:3);%提取前三个主成分得分值 sc%输出前三个主成分得分值 e=eigenvalue(1:3)';%提取前三个特征根并转置 M=e(ones(m,1),:).^0.5;%输出前三个特征根并转置 compmat=p3.*M;%利用特征根构造变换矩阵 per=100*eigenvalue/sum(eigenvalue);%求出成分载荷矩阵的前三列 per %求出各主成分的贡献率 cumsum(per);%列出各主成分的累积贡献率 figure(1) pareto(per);%将贡献率绘成直方图 t2 figure(2) %输出各省与平局距离 plot(eigenvalue,'r+');%绘制方差贡献散点图 hold on %保持图形 plot(eigenvalue,'g-');%绘制方差贡献山麓图

用MATLAB实现结构可靠度计算.

用MATLAB实现结构可靠度计算 口徐华…朝泽刚‘u刘勇‘21 。 (【l】中国地质大学(武汉工程学院湖北?武汉430074; 12】河海大学土木工程学院江苏?南京210098 摘要:Matlab提供了各种矩阵的运算和操作,其中包含结构可靠度计算中常用的各种数值计算方法工具箱,本文从基本原理和相关算例分析两方面,阐述利用Matlab,编制了计算结构可靠度Matlab程.序,使得Matlab-语言在可靠度计算中得到应用。 关键词:结构可靠度Matlab软件最优化法 中图分类号:TP39文献标识码:A文章编号:1007-3973(200902-095-Ol 1结构可靠度的计算方法 当川概率描述结构的可靠性时,计算结构可靠度就是计算结构在规定时问内、规定条件F结构能够完成预定功能的概率。 从简单到复杂或精确稃度的不同,先后提出的可靠度计算方法有一次二阶矩方法、二次二阶矩方法、蒙特卡洛方法以及其他方法。一次■阶矩方法又分为。I-心点法和验算点法,其中验算点法足H前可靠度分析最常川的方法。 2最优化方法计算可靠度指标数学模型 由结构111n个任意分布的独立随机变量一,x:…以表示的结构极限状态方程为:Z=g(■.托…t=0,采用R-F将非正念变量当罱正态化,得到等效正态分布的均值o:和标准差虹及可靠度指标B,由可靠度指标B的几何意义知。o;辟

开始时验算点未知,把6看成极限状态曲面上点P(■,爿:---37,的函数,通过优化求解,找到B最小值。求解可靠皮指标aJ以归结为以下约束优化模型: rain睁喜t华,2 s.,.Z=g(工i,x2’,…,工:=0 如极限状态方栉巾某个变最(X。可用其他变量表示,则上述模型jfIJ‘转化为无约束优化模型: 。。B!:手f生丛r+阻:坚:坠:盐尘}二剐 t∞oY?’【叫,J 3用MATLAB实现结构可靠度计算 3.1Matlab简介 Matlab是++种功能强、效率高、便.丁.进行科学和工程计算的交互式软件包,汇集了人量数学、统计、科学和工程所需的函数,MATI.AB具有编程简甲直观、用户界mf友善、开放性强等特点。将MATLAB用于蒙特卡罗法的一个显著优点是它拥有功能强大的随机数发生器指令。 3.2算例 3.2.I例:已知非线形极限状态方程z=g(t r'H=567f r-0.5H2=0’f、r服从正态分布。IIf=0.6,o r=0.0786;la|_ 2.18,o r_0.0654;H服从对数正态分布。u H= 3218,O。 =0.984。f、r、H相互独立,求可靠度指标B及验算点(,,r’,H‘。 解:先将H当量正念化:h=ln H服从正态分布,且 ,‘-““了:等专虿’=,。49?口二-、『五ir面_。。3

BP神经网络matlab源程序代码

BP神经网络matlab源程序代码） %******************************% 学习程序 %******************************% %======原始数据输入======== p=[2845 2833 4488;2833 4488 4554;4488 4554 2928;4554 2928 3497;2928 3497 2261;... 3497 2261 6921;2261 6921 1391;6921 1391 3580;1391 3580 4451;3580 4451 2636;... 4451 2636 3471;2636 3471 3854;3471 3854 3556;3854 3556 2659;3556 2659 4335;... 2659 4335 2882;4335 2882 4084;4335 2882 1999;2882 1999 2889;1999 2889 2175;... 2889 2175 2510;2175 2510 3409;2510 3409 3729;3409 3729 3489;3729 3489 3172;... 3489 3172 4568;3172 4568 4015;]'; %===========期望输出======= t=[4554 2928 3497 2261 6921 1391 3580 4451 2636 3471 3854 3556 2659 ... 4335 2882 4084 1999 2889 2175 2510 3409 3729 3489 3172 4568 4015 ... 3666]; ptest=[2845 2833 4488;2833 4488 4554;4488 4554 2928;4554 2928 3497;2928 3497 2261;... 3497 2261 6921;2261 6921 1391;6921 1391 3580;1391 3580 4451;3580 4451 2636;... 4451 2636 3471;2636 3471 3854;3471 3854 3556;3854 3556 2659;3556 2659 4335;... 2659 4335 2882;4335 2882 4084;4335 2882 1999;2882 1999 2889;1999 2889 2175;... 2889 2175 2510;2175 2510 3409;2510 3409 3729;3409 3729 3489;3729 3489 3172;... 3489 3172 4568;3172 4568 4015;4568 4015 3666]'; [pn,minp,maxp,tn,mint,maxt]=premnmx(p,t); %将数据归一化 NodeNum1 =20; % 隐层第一层节点数 NodeNum2=40; % 隐层第二层节点数 TypeNum = 1; % 输出维数 TF1 = 'tansig';

计算方法上机实验报告-MATLAB

《计算方法》实验报告 指导教师： 学院： 班级： 团队成员：

一、题目 例2.7应用Newton 迭代法求方程210x x --=在1x =附近的数值解 k x ，并使其满足8110k k x x ---< 原理： 在方程()0f x =解的隔离区间[],a b 上选取合适的迭代初值0x ，过曲线()y f x =的点()() 00x f x ，引切线 ()()()1000:'l y f x f x x x =+- 其与x 轴相交于点：()() 0100 'f x x x f x =-，进一步,过曲线()y f x =的 点()()11x f x , 引切线 ()()()2111: 'l y f x f x x x =+- 其与x 轴相交于点：() () 1211 'f x x x f x =- 如此循环往复，可得一列逼近方程()0f x =精确解*x 的点 01k x x x ,,,,,其一般表达式为： ()() 111 'k k k k f x x x f x ---=- 该公式所表述的求解方法称为Newton 迭代法或切线法。

程序： function y=f(x)%定义原函数 y=x^3-x-1; end function y1=f1(x0)%求导函数在x0点的值 syms x; t=diff(f(x),x); y1=subs(t,x,x0); end function newton_iteration(x0,tol)%输入初始迭代点x0及精度tol x1=x0-f(x0)/f1(x0);k=1;%调用f函数和f1函数 while abs(x1-x0)>=tol x0=x1;x1=x0-f(x0)/f1(x0);k=k+1; end fprintf('满足精度要求的数值为x(%d)=%1.16g\n',k,x1); fprintf('迭代次数为k=%d\n',k); end 结果：

计算方法及其MATLAB实现第二章作业

作者：夏云木子 1、 >> syms re(x) re(y) re(z) >> input('计算相对误差：'),re(x)=10/1991,re(y)=0.0001/1.991,re(y)=0.0000001/0.0001991 所以可知re(y)最小，即y精度最高 2、 >> format short,A=sqrt(2) >> format short e,B=sqrt(2) >> format short g,C=sqrt(2)

>> format long,D=sqrt(2) >> format long e,E=sqrt(2) >> format long g,F=sqrt(2) >> format bank,H=sqrt(2) >> format hex,I=sqrt(2) >> format +,J=sqrt(2) >> format,K=sqrt(2)

3、 >> syms A >> A=[sqrt(3) exp(7);sin(5) log(4)];vpa(pi*A,6) 4、1/6251-1/6252=1/6251*6252 5、（1）1/（1+3x）-(1-x)/(1+x)=x*(3*x-1)/[(1+3*x)*(1+x)] (2) sqrt(x+1/x)-sqrt(x-1/x)=2/x/[sqrt(x-1/x)+sqrt(x+1/x)] (3) log10(x1)-log(x2)=log10(x1/x2) (4) [1-cos(2*x)]/x =x^2/factorial(2)-x^4/factorial(4)+x^6/factorial(6)-…