8、数值分析与matlab

1、今天要拷matlab代码了，而且是很恶心的算法，估计也没几个人能看得懂，就连我自己都看不懂。

　　我也不知道这样做的意义何在，可能只是证明我在这世上曾经学过那么那么难的东西吧

　　首先是一个matlab版的快速排序，同学们应该都看得懂吧。

function f=quicksort(x,left,right)
if left<right
    [i,x]=Division(x,left,right);
    x=quicksort(x,left,i-1);
    x=quicksort(x,i+1,right);
end
f=x;

function [i,x]=Division(x,left,right)
base=x(left,2);
while left<right
    while left<right&x(right,2)>=base
        right=right-1;
    end
    c=x(left,2);d=x(right,2);
    c1=x(left,1);d1=x(right,1);
    c2=x(left,3);d2=x(right,3);
    c3=x(left,4);d3=x(right,4);
    c5=x(left,5);d5=x(right,5);
    x(left,2)=d;x(right,2)=c;
    x(left,1)=d1;x(right,1)=c1;
    x(left,3)=d2;x(right,3)=c2;
    x(left,4)=d3;x(right,4)=c3;
    x(left,5)=d5;x(right,5)=c5;
%     x(left,1)=x(right,1);
    while left<right&x(left,2)<=base
        left=left+1;
    end
    c=x(left,2);d=x(right,2);
    c1=x(left,1);d1=x(right,1);
    c2=x(left,3);d2=x(right,3);
    c3=x(left,4);d3=x(right,4);
    c5=x(left,5);d5=x(right,5);
    x(left,2)=d;x(right,2)=c;
    x(left,1)=d1;x(right,1)=c1;
    x(left,3)=d2;x(right,3)=c2;
    x(left,4)=d3;x(right,4)=c3;
    x(left,5)=d5;x(right,5)=c5;
%     x(right,1)=x(left,1);
end
i=left;

　　以上大概的意思就是根据向量中第二列的值，将向量的其他列进行快速排序

2、下边这个应该是二分法求函数零点的程序吧

function [xstar,index,it]=bisect(fun,a,b,ep)
if nargin<4 ep=1e-5;end
fa=feval(fun,a);
fb=feval(fun,b);
if fa*fb>0
    xstar=[fa,fb];index=0;it=0;
    return
end
k=0;
while abs(b-a)/2>ep
    x=(a+b)/2;fx=feval(fun,x);
    if fx*fa<0
        b=x;fb=fx;
    else
        a=x;fa=fx;
    end
    k=k+1;
end
xstar=(a+b)/2;index=1;it=k;

3、逆天的chi2plot,也就是传说中的正态概率图，属于数据分析部分

function chi2plot(X)
dd=[];
p=[];
[M,N]=size(X);
MEAN=mean(X);
SS_1=inv(cov(X));
for byk=1:M;
    DD=(X(byk,:)-MEAN)*SS_1*(X(byk,:)-MEAN)';
    dd=[dd,DD];
    pp=(byk-0.5)/M;
    p=[p,pp];
end
dd=sort(dd)'
xx=chi2inv(p,N)'
plot(xx,dd,'+'),lsline
xlabel('chi2quantitle')
ylabel('Sample generalized diatance')
title('chi2plot')

4、改进的欧拉公式

function [x,y]=Euler_correct(fun,a,b,n,y0)
%改进的Euler公式，其中
%fun为一阶微分方程的函数
%a,b为求解区间的左右端点
%n为等分区间；
%y0为初始条件
x=zeros(1,n+1);y=zeros(1,n+1);
h=(b-a)/n;
x(1)=a;y(1)=y0;
for k=1:n
    x(k+1)=x(k)+h;
    y0=y(k)+h*feval(fun,x(k),y(k));
    y(k+1)=y(k)+h/2*(feval(fun,x(k),y(k))+feval(fun,x(k+1),y0));
end