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);endf=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);endi=left;

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

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

function [xstar,index,it]=bisect(fun,a,b,ep)if nargin<4 ep=1e-5;endfa=feval(fun,a);fb=feval(fun,b);if fa*fb>0    xstar=[fa,fb];index=0;it=0;    return endk=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;endxstar=(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];enddd=sort(dd)'xx=chi2inv(p,N)'plot(xx,dd,'+'),lslinexlabel('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