Matlab数值微分

 1 %LU分解算法
 2 %输入：矩阵A和向量b
 3 %输出：det和自变量的值
 4 function [det,x,y,l,u]=LU(A,b)
 5 [m,n]=size(A);
 6 %输入的A与b的行数与列数的限制
 7 if m~=n
 8     disp(‘输入错误，系数矩证阵只能是方阵‘);
 9 end
10 if n~=length (b)
11     disp(‘输入错误，常数项的个数应与方程的个数相同‘);
12 end
13 for k=1:n-1%主行循环
14     for i=k+1:n
15         A(i,k)=A(i,k)/A(k,k);
16         for j=k+1:n
17             A(i,j)=A(i,j)-A(i,k)*A(k,j);
18         end
19     end
20 end
21 l=eye(n);
22 u=zeros(n,n);
23 for i=1:n
24     for j=1:n
25         if j<i
26             l(i,j)=A(i,j);
27         else
28             u(i,j)=A(i,j);
29         end
30     end
31 end
32 
33 y(1)=b(1);
34 for k=2:n
35     s=0;
36     for j=1:k-1
37         s=s+l(k,j)*y(j);
38     end
39     y(k)=b(k)-s;
40 end
41 x(n)=y(n)/u(n,n);
42 for k=n-1:-1:1
43     s=0;
44     for j=k+1:n
45         s=s+u(k,j)*x(j);
46     end
47     x(k)=(y(k)-s)/u(k,k);
48 end
49 det=1;
50 for i=1:n
51     det=det*u(i,i);
52 end
53 end

 1 %列主元消去算法，求自变量和系数行列式的值
 2 %输入：系数矩阵和因变量，也就是Ax=b,中的A和b
 3 %输出：自变量x和系数行列式det
 4 function [z,det]=liezhu(A,b)
 5 %行列式det初始值为1
 6 det=1;
 7 [m,n]=size(A);
 8 %输入的A与b的行数与列数的限制
 9 if m~=n
10     disp(‘输入错误，系数矩证阵只能是方阵‘);
11 end
12 if n~=length (b)
13     disp(‘输入错误，常数项的个数应与方程的个数相同‘);
14 end
15 %对于k=1,2,...,n-1，对A扫描
16 for k=1:n-1
17     %按列选主元
18     p=A(k,k);%记忆当前值
19     I=k;
20     for i=k:n
21         if abs(A(i,k)) > abs(p)%按行扫描该列的最大值
22         end
23     end
24     if I~=k%换行使当前主元绝对值最大
25         for j=k:n
26             w=A(k,j);%记忆当前主元所在行
27             A(k,j)=A(I,j);%换行
28             A(I,j)=w;%换行
29         end
30         %并且把b同时换行
31             u=b(k);
32             b(k)=b(I);
33             b(I)=u;
34     end
35     %如果得到的最大绝对值是0，则结束程序
36     if A(i,k)==0
37         det=0;
38     end
39     %消元计算
40     for i= k+1:n
41     %对于i=k+1,...,n
42         A(i,k)=A(i,k)/A(k,k);
43         b(i)=b(i)-A(i,k)*b(k);
44         for j=k+1:n
45         %对于j=k+1,...,n
46             A(i,j)=A(i,j)-A(i,k)*A(k,j);
47         end
48     end
49     %算行列式
50     det=A(k,k)*det;
51 end
52 %如果最后一个主元等于0，则停止计算。并使行列式等于0
53 if A(n,n)==0
54     det=0;
55 end
56 %回代求解
57 b(n)=b(n)/A(n,n);
58 %对于i=n-1,...,2,1
59 for i=(n-1):-1:1
60     w=0;
61     %得到一个累加值
62     for j=(i+1):n
63         w=w+A(i,j)*b(j);
64     end
65     %求得bi
66     b(i)=(b(i)-w)/A(i,i);
67 end
68 %行列式的值
69 det=det*A(n,n);
70 z=b; 
71 end