数值分析综述报告

淮阴工学院

《数值分析》考试

──基于Matlab的方法综合应用报告

班级：金融1121 姓名：姚婷婷

学号： 1124104129

成绩：

数理学院

2014年6月7日

《数值分析》课程综述报告

前言：

数值分析也称计算方法，它与计算工具的发展密切相关。数值分析是一门为科学计算提供必需的理论基础和有效、实用方法的数学课程，它的任务是研究求解各类数学问题的数值方法和有关的理论。 正文：

第一章 近似计算与误差分析

1、产生误差的原因：①模型误差；②观测误差；③截断误差；④舍入误差。

2、四则运算的误差： ①加减法运算

()()()****x y x y δδδ±=+

②乘法运算

()()()

******

*********xy x y xy xy xy x y x y y y x x x y x y y x δδ

δ-=-+-≤-+-?=+

③ 除法运算：

()()()

()

*****

*******

*

*

**

*

*

*

*

*

**

x x xy x y y y yy xy x y x y x y

yy x x y

y y x yy x y

y x

x y y δδ

δ--=-+-=-+-=

+??

?≈

???

3、科学表示法、有效数字、近似值的精度 任何一个实数都可以表示成如下的形式：

其中：是正整数，是整数，

如果是数的近似值

并且

则称该近似值具有位有效数字（significant digit ）。 此时，该近似值的相对误差为

另一方面，若已知

()()

*111

1021n

r x a δ-≤

+

那么，

()()***1112110.10

211

102

r m n n m n x x x x a a a a δ----≤?=+≤

即：*x 至少有n 位有效数字。 例： 3.141592653589793...π= 取其近似值如下： x*=3.14 x *

=3.14159 x*=3.1415 x*=3.141

*

*213

100.314

110.0016...0.00510102

2

x x π--=?-=<=?=?

**516

100.314159

110.0000026...0.00000510102

2

x x π--=?-=<=?=?

*

*314

100.31415

110.000092...0.0001101022

x x π--=?-=<

*

*213

100.3141

110.00059...0.001101022

x x π--=?-=<

第二章 线性方程组

在科学计算中，问题的本身就是求解线性方程组，许多问题的求解需要最后也归结为线性方程组的求解，所以线性方程组的求解是科学计算中最常见的问题。

对于线性方程组的求解一般有两种方法：（1) 直接法：高斯消去法；（2) 间接法：各种迭代法。

(1) 高斯消去法：

①求解思路：先消元，即按一定的规律逐步消去未知量，将方程组Ax B =化为等价的上（或下）三角形方程组；然后进行回代，即由上三角形方程组逐个求出；

②高斯（列、全）主元素消去法,及在消元的每一步选取（列）主元素——列中绝对值最大的元取做主元素，计算步骤：⑴消元过程：按列选主元、行交换、消元计算；⑵回代过程； ③高斯列主元素消去法的MATLAB 实现：。

第三章 解线性方程组的迭代法

通常逆矩阵不易求得，特别是对于大型的线性方程组，需要用迭代法求解。

用迭代法求解线性方程组，要把线性方程组写成等价的形式，右边写为迭代格式，如：

k

b x k A k x b x k A b kx Ax kx b Mx b x E A b x Ax x b Mx x b Ax E E n n ++=

∴++=++=?+=++=++=?+=?-=)(1)()(

2、关于迭代法收敛性的两个重要结论： ①充分必要条件是：矩阵的谱半径

()1

M ρ<；

②充分条件是：矩阵M 的某个算子范数

1

M <。

3、线性方程组的迭代法主要有Jacobian 迭代法，Gauss-Seidel 迭代法。 ①Jacobian 迭代法:

()()(

)

()()

()

11111

k k Ax b D L U x

b

A D L U

Dx L U x b x D L U x D b

M D L U f D b

+----=?

?--=?

=--??=++?=++?=+??=?

②Gauss-Seidel 迭代法:

()(

)

()(

)

()()()11

11

1

k k Ax b D L x Ux b A D L U x D L Ux D L b

M D L U f D L b

--+--=?

?-=+?

=--?=-+-=-?

=- (3.7)

③Jacobian 迭代法与G-S 迭代法比较：

()()()()()()()

()()111111*********

1121112100

00

00000000k k k k k k n nn n n k n k n n k n x x a x x D a a x x x a a x D D b

a x ++++-++----?????? ? ???- ? ???

= ? ??? ? ??? ? ?--??????

??

--?? ?

?? ???++ ?-?? ??? ??

??

?

（3.8）

式（3.7）和 (3.8) 表明：Gauss-Seidel 迭代法在计算第1k

+次迭代的第i 个分量()1k i

x +时，

及时地利用了在此步迭代中得到的新的迭代值：()1k j

x +，

1,2,

,1j i =-，由于第1k +步的迭

代值通常比第k 步的迭代值更接近方程组的精确解，所以，在Jacobian 迭代法和GS 迭代法都收敛的情况下，Gauss-Seidel 迭代法的收敛速度比Jacobian 迭代法的收敛速度高。

例题：用MATLAB 函数normrdn 生成5阶矩阵M 和向量b 分别构造线性方程组Ax b =的Jacobi 迭代格式和G-S 迭代格式，并判断收敛性。

Jacobian 迭代法和GS 迭代法程序如下：

clc;

clear all;

%1、生成M 和b M=normrnd(1,2,5) b=normrnd(1,2,5,1) %Jacobian 迭代法 M1=D\(L+U) f1=D\b

rho=max(abs(eig(M1)));

R=1e-08; %设定的一个收敛标准 switch sign(1-rho) case -1

disp('the Jocobian method is not applicable') otherwise

x(:,1)=normrnd(0,9,5,1); k=1 while k<=50*5

x(:,k+1)=M1*x(:,k)+f1; if norm(x(:,k+1)-x(:,k))>=R k=k+1; else

X=x(:,k+1);

disp('Jacobian 迭代法迭代次数为：') IterN=k %Jacobian 迭代法迭代次数 break end end end

%Causs-Seidel 迭代法 M2=(D-L)\U f2=(D-L)\b

rho=max(abs(eig(M2))); R=1e-08;

switch sign(1-rho) case -1

disp('the auss-seidel method is not applicable')

otherwise

x(:,1)=normrnd(0,9,5,1); k=1 while k<=50*5

x(:,k+1)=M2*x(:,k)+f2; if norm(x(:,k+1)-x(:,k))>=R k=k+1; else

X=x(:,k+1)

disp('Causs-Seidel 迭代法迭代次数为：') IterN=k break end end end

第四章 非线性代数方程（组）的数值解法： 一、二分法：

首先要确定适当的包含根的区间，这可以依据闭区间上连续函数的介值定理来确定，例如该方程：

()320.80.750

f x x x x =--+=

x

f (x )

对于该方程

()120.80.750

f -=-++<

()00.750

f =>

所以该方程至少有一个实根位于区间，图像表明该区间中只含有一个实根；

用*

x 表示方程()0

f x =在区间

[],a b 上的精确解，对于给定的精度要求0ε>，取区间[]

,a b 的中点

2

a b x +=

，并按下式进行判断：

()()()()()***00[,]0[,]

f x x x f x f a x a x f x f b x x b ?=?=??

?

?

（1）

以

()()0

f x f a <为例，

① 如果

2

b a

ε-≥,没有达到精度要求，令x

b =，并重复（1）的迭代过程；

②如果

2

b a

ε-<，那么，必有

[]

*,,x x x a x ε-

，

因为

[]*,x a x ∈。即区间[]

,a x 内的任何一点都可以作为方程()0f x =的近根，特别地，可取x 做为近似解。

二、牛顿迭代法：

任取初始值)(],,[0x f y b a x =∈上过点（)(,0

0x x f ）的切线方程为：))(()(000x x x x f f y -'+=与

x 轴交于点x 1，

)

()

(0001x x x x f f '-

=，过点))(,(1

1

x x f 的切线方程为))(()(1

1

1

x x x x f f y -'+=与x 轴交于点

x 2

，

)()

(1112x x x x f f '-

= ，……… ，如此下去得牛顿迭代公式：

)

()

(1

x x x

x

n n n

n f f '-

=

+

例题：考虑如下三阶非线性方程组：

()()()222

22222

=2

220

4a x y z a x y a a x a y a z b ì?++??

+=í??-+-+-=??

其中

取适当的迭代初值()000,y ,T

x z ，用Newton 迭代法求方程组的数值解

程序：

%Newton 迭代法求解 x=sym('x','clear'); y=sym('y','clear'); syms z ;

F=[x^2+y^2+a^2*z^2/2-a^2; x+y-a;

(2*x-a)^2+(2*y-a)^2+a^2*(z-b)/4]; Fx=diff(F,x,1); Fy=diff(F,y,1); Fz=diff(F,z,1); DF=[Fx Fy Fz];

F=@(x,y,z)[x^2+y^2+a^2*z^2/2-a^2;x+y-a; (2*x-a)^2+(2*y-a)^2+a^2*(z-b)/4];

%Newton 迭代法求解过程 Fr=10^-10; Er=10^-10; x0=[-a/4;-a;-a/2]; x0=[-a/4;-a;a/2]; k=1; X(:,1)=x0;

while norm(F(X(1,k),X(2,k),X(3,k))-[0;0;0],2)>=Fr tic;

f(:,k)=F(X(1,k),X(2,k),X(3,k)); J=subs(DF,'x',X(1,k)); J=subs(J,'y',X(2,k));

J=subs(J,'z',X(3,k)); X(:,k+1)=X(:,k)-J\f(:,k); t(k)=toc;

if norm((X(:,k+1)-X(:,k)),2)

disp('Newton 迭代法结果为：'); disp(X(:,end)); 运行结果：

Newton 迭代法结果为： 3.4291

0.4658

0.6535

第五章 插值

一、插值：插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。

二、常用差值法：拉格朗日（Lagrange ）插值法、牛顿（Newton ）插值法 1、拉格朗日（Lagrange ）插值法： 拉格朗日插值多项式

()()0

n

k k

k L x L x y ==∑ ⑴简单的证明

因为拉格朗日插值多项式的基函数有如下的性质：

()()()01

,n

j k

k j j k

j

j k

x x x x L x x x j k

x x =≠-=?==?

=≠-?∏

0,1,2,

,k n =

所以拉格朗日插值多项式

()()()0n

k j k j k k k k

j L x L x y L x y y ====∑ 0,1,2,,k n =

满足插值的条件。 插值多项式

()()0n

k k

k L x L x y ==∑

⑶拉格朗日插值法的不足

在实际问题中，观测的数据可能会不断增加，如果用拉格朗日插值公式构造插值多项式，那么，每当增加数据就要重新计算多项式的系数，由此增加许多不必要的计算工作量。

2、（三次）样条（Spline ）插值 ⑴插值条件要求

要求所求的插值多项式()S x （三次样条函数） ① 在每个区间[

]1,k k x x -，0,1,2,,k n =，是次数不超过三次的多项式；

②()k k S x y =，0,1,2,

,k n =；

③ ()S x 在区间..上具有二阶连续导数。

例题：在夏季，较大湖泊的水体按深度被跃变层分为上部的变温层和下部的均温层。水体的分层化对环境工程中污染问题的研究具有重要的意义，例如，有机物的分解将导致被跃变层隔离的底部水体中氧急剧减少。按温度随水深的变化曲线()x T T =，跃变层位于水深处：

()00max x bottom

d x T x dx ≤≤??

=????

现有美国普拉特湖（Platte Lake ） 的一组数据：

1. 试利用Lagrange 差值方法求温度随水深（近似）变化函数

()

l T T x =表达式；

2. 试利用三次样条差值方法（应用Matlab 函数csape ）求温度随水深（近似）变化函数

()s T T x =表达式；

3. 画出插值函数()l T T x =和()s T T x =曲线，并与原始插值数据图像作比较 程序代码：

函数文件程序：

function Ln = my_Fun(x, XI, YI)

if isa(x, 'sym') == 1;

n = length(XI) - 1;

Ln = 0;

Pn = sym(ones(n + 1, 1));

for k = 1 : n + 1

for i = 1 : n + 1

if i ~= k

Pn(k) = Pn(k)*(x - XI(i))/(XI(k) - XI(i));

else

Pn(k) = Pn(k);

end

end

Ln = Ln + Pn(k)*YI(k);

end

else

n = length(XI) - 1;

L = ones(n + 1, length(x));

Ln = zeros(size(x));

for k = 1 : n + 1

for i = 1 : n + 1

if i ~= k

L(k,:) = L(k,:).*(x - XI(i))./(XI(k) - XI(i));

else

L(k,:) = L(k,:);

end

end

Ln = Ln + L(k, :).*YI(k);

end

End

主文件程序：

clc;

clear all;

close all;

%问题1.用Lagrange差值方法求温度随水深（近似）变化图像

t = sym('t', 'real');

x=[0 2.3 4.9 9.1 13.7 18.3 22.9 27.2];

T=[22.8 22.8 22.8 20.6 13.9 11.7 11.1 11.1];

X=linspace(0,27.2,275);

Ln=my_Fun(X,x,T);

figure(1)

plot(x,T,'r*','markersize',15)

xlabel('深度x');

ylabel('温度T');

title('原始的散点图');

pause(3)

hold on

plot(X,Ln,'b-','linewidth',3);

xlabel('深度x');

ylabel('温度T');

title('Lagrange差值图');

%求出温度随水深（近似）变化函数表达式

I_poly = my_Fun(t, x, T);

I_poly = simple(I_poly)

I_poly = sym2poly(I_poly)

%2.用三次样条差值方法求温度随水深（近似）变化函数表达式

figure(2)

hold on

y=interp1(x,T,X,'spline');

plot(X,y,'b-','linewidth',3); %用三次样条插值方法画出图像

xlabel('深度x');

ylabel('温度T');

title('三次样条法差值图');

PP=csape(x,T,[2,2],[0,0]);

Coefs = PP.coefs

%3.两种差值函数图象比较

figure(3)

plot(X,Ln,'-','color',[1, 0, 0],'LineWidth', 3)

xlabel('深度x');

ylabel('温度T');

title('Lagrange差值图');

hold on

pause(3)

fnplt(PP,'b',3,[0,28]) %函数作图

xlabel('深度x');

ylabel('温度T');

title('三次样条法差值图');

pause(3)

plot(x,T,'--','linewidth',3,'markersize',10)

运行结果:

I_poly =

(11156832321404503197312500000*t^7 - 942692507666801544174224450000*t^6 + 29858157926737187739430003070000*t^5 - 433725749060135044183186809635000*t^4 + 2850751947133625442244870404829250*t^3 -

8285950923799759686915051372346805*t^2 + 8477820147617239232740036583612328*t + 512657632901869980094703083834211328)/22484983899204823688364170343605760

I_poly =

Columns 1 through 6

0.0000 -0.0000 0.0013 -0.0193 0.1268 -0.3685 Columns 7 through 8 0.3770 22.8000 Coefs =

0.0022 -0.0000 -0.0115 22.8000 -0.0092 0.0150 0.0230 22.8000 -0.0114 -0.0565 -0.0850 22.8000 0.0297 -0.2004 -1.1640 20.6000 -0.0153 0.2099 -1.1200 13.9000 0.0017 -0.0014 -0.1605 11.7000 -0.0017 0.0223 -0.0640 11.1000

图示如下：

深度x

温度T

原始的散点图

深度x

温度T

三次样条法差值图

深度x

温度T

Lagrange 差值图

深度x

温度T

三次样条法差值对比图

第六章 最小二乘拟合与最佳逼近 一、最小二乘拟合

①加权最小二乘法逼近准则：

()min ,P d f P ∈P

=min

P ∈P

()()1

20N

i i i i f x P x ρ=????-??????

∑

②最小二乘逼近多项式

()

*n P x 必须满足如下必要条件：

1

220T T n a

a f a a

a

????Φ?? ??

??Φ

?

?Φ

??==-+= ??

? ??Φ ? ????

即满足法方程组:

()1

T T T T a f

a f

??????-=?=

试利用最小二乘法求满足

()()()()1

22*

2

00.5min n N n k n k P x k f P f x P x m

=??-=-≤????

∑ 即误差不超过0.5m 的最低次数的拟合多项式()*n P x ,写出该多项式的表达式;画出数据散点图和该多项式曲线.

程序：

clc; clear all ; close all ;

x=linspace(0,12,12);

x=[0 200 1000 2100 3500 5000 6800 7500 9000 11200 12000]; h=[1.64 1.58 1.68 1.84 1.58 0.86 0.39 0.31 0.39 0.77 0.86]; plot(x,h,'*','markersize',8) %画出给出数据的散点图 xlabel('水平距离x'); ylabel('高度h');

title('岩层断面水平距离x 和高度h 的散点图'); figure(1)

%1.最小二乘法拟合 n=3;

[P, S ] = polyfit(x, h, n) Pm = polyval(P, x)

R = sqrt(sum((h - Pm).^2)) %误差

t = linspace(x(1), x(end), 12); Poly = polyval(P, t); figure(2)

plot(x, h, 'ro', t, Poly, 'LineWidth', 3, 'markersize', 8) set(gca,'FontSize',18)

legend('The Data', 'The Fitting Curve', 1)

title('Curve Fitting by Least Square Approximation', 'fontsize', 18)

水平距离x

高度h

岩层断面水平距离x 和高度h 的散点图

Curve Fitting by Least Square Approximation

第七章 微积分的数值方法 一、数值微分

如果给定函数 的关系式比较复杂或者 未知，而仅仅知道在1n +个相异点k

x ，0,1,,k

n

=处的函数值，那么，我们可以利用函数的插值多项式的导数作为函数导数的近似值

()()()()

0n

n k k k f x L x l x f x =≈=∑

因而有

()()n

f x L x ''≈

这里需要说明一点的是，尽管和的函数值可能相差不多，但是它们的导数有可能相差很大。 二、数值积分

数值积分方法最大的优点是将复杂的函数积分转化为相对简单的加、减、乘、除运算。

对于给定的函数，如果 的函数关系式比较复杂或.未知，而仅仅知道该函数在1n +个点处的函数值k y ，

那么可以利用函数的插值多项式的积分作为函数在

[,]a b 上的定积分的近似值。

1、Newton-Cotes 求积公式 根据Lagrange 插值多项式

()()()

0n

n k k k L x l x f x ==∑

有 ()()()()0b

a

b

n a

n

b

k k a k I f x dx

L x dx

l x dx f x ==≈??=??????∑?

令

()b

k k a

A l x dx

=? 0k n ≤≤

得Newton-Cotes 求积公式()()()()()

00

b

b

n a

a

b

n

k k k a

n

k k k f x dx L x dx

f x l x dx

A f x ==≈==?

?∑?∑

特别地，当取插值节点为,0,1,,k x a kh k n

b a h n

=+=-=

时有

①两点公式（梯形公式）：

()()()

()()0

2b

b

n a

a

n

k k k f x dx L x dx

A f x b a

f a f b =≈=-??=

+?

??

?∑

②三点公式（Simpson 公式）：()()()()462b

b

n a

a

f x dx L x dx

b a a b f a f f b ≈?

?-+??

=

++ ???????

?

?

0,12内所围⑴试利用复化的梯形求积法求该组数据所在曲线与基准线（x轴）在范围[]

成图形面积.画出数据散点图和图形的示意图.

0,9内所⑵试利用复化的Simpson求积法求该组数据所在曲线与基准线（x轴）在范围[]

围成图形面积.画出数据散点图和图形的示意图.

解答如下：

程序：

% 题目给出的数据

XI = [0 0.20 1.00 2.10 3.50 5.00 6.80 7.50 9.00 11.2 12.00];

YI = -[1.64 1.58 1.68 1.84 1.58 0.86 0.39 0.31 0.39 0.77 0.86];

n = length(XI) - 1; %此时n=10

M = 60;

x = linspace(XI(1), XI(end), M +1); %设置线性空间

% 问题1.1

T1 = 0; %初值

%将图像分成若干个小区间，在每个区间内，求小梯形的面积，并累加起来，就是题目所要求的面积

for k = 2 : n + 1

T1 = T1 + (YI(k) + YI(k-1))*(XI(k) - XI(k-1))/2;

end

disp('复化的梯形求积法求得的图形面积：')

T_Trape = T1 %累加后的图形面积

figure(1) %画出函数图像

set(gca,'FontSize', 20)

patch([XI(1), XI, XI(end)], [0, YI, 0], [0.5 1 0.5]) %颜色的填充

hold on

plot(XI, YI, 'o-', XI, 0*XI, 'k', 'LineWidth', 3, 'markersize', 8)

title('复化的梯形求积法求围成图形面积')

xlabel('水平距离x')

ylabel('深度y')

运行结果：

复化的梯形求积法求得的图形面积：

T_Trape =

-11.6090

图像：

复化的梯形求积法求围成图形面积

水平距离x

深度y

程序：

函数M 文件：

function Ln = Lagrange_Fun_01(x, XI, YI) if isa(x, 'sym') == 1; n = length(XI) - 1; Ln = 0;

Pn = sym(ones(n + 1, 1)); for k = 1 : n + 1 for i = 1 : n + 1 if i ~= k

Pn(k) = Pn(k)*(x - XI(i))/(XI(k) - XI(i)); else

Pn(k) = Pn(k); end end

Ln = Ln + Pn(k)*YI(k); end else

n = length(XI) - 1;

L = ones(n + 1, length(x)); Ln = zeros(size(x)); for k = 1 : n + 1 for i = 1 : n + 1 if i ~= k

L(k,:) = L(k,:).*(x - XI(i))./(XI(k) - XI(i)); else

L(k,:) = L(k,:); end end

Ln = Ln + L(k, :).*YI(k); end end 主程序：

% 问题1.2 和 1.3

%问题3是问题2的延伸 S2 = 0;

%Simpson 复化求积法，应用同样的面积累加，得到最终的图形面积 for k = 1: 2: n - 1

S2 = S2 + (XI(k+2) - XI(k))*(YI(k) + 4*YI(k+1) + YI(k+2))/6; X(k, :) = linspace(XI(k), XI(k+2), 20);

Y(k, :) = Lagrange_Fun_01(X(k,:), XI(k : k+2), YI(k : k+2)); end

disp('Simpson 复化求积法求得的图形面积：') S_Simpson = S2 %Simpson 复化求积法求得的面积 figure(2)

set(gca,'FontSize', 20) hold on

for k = 1: 2 : n-1

patch([XI(k), X(k,:), XI(k +2)], [0, Y(k, :), 0], [0.5 1 0.5]) plot(X(k, :), Y(k, :), 'LineWidth', 3) end

plot(XI, YI, 'o', XI, 0*XI, 'k', 'LineWidth', 3, 'markersize', 8) title('复化的Simpson 求积法求围成图形面积') xlabel('水平距离x') ylabel('深度y') 运行结果：

Simpson 复化求积法求得的图形面积： S_Simpson = -11.9128 图像：

复化的Simpson 求积法求围成图形面积水平距离x

深度y

第八章 微分方程（组）初值问题数值方法

引言：微分方程（组）是科学研究和工程应用中最常用的数学模型之一。但是基绝大多数的常微分方程和常微分方程组得不到解析解和实际应用中往往只需要知道常微分方程（组）的解在某些点处的函数值这两点原因，绝大多数在实际中遇到的常微分方程和常微分方程组得

不到“解析解”，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。

一、 Euler 方法——最简单的数值解法 1、Euler 法 假设要求在点

0k t t kh

=+,，1,2,,k n =

处初值问题（4）的解的近似值。首先对式（4）的两端积分，得

()()()()11

1,k k k k t t k k t t dy dt y t y t f t y t dt

dt +++=-=??

对于该式的右边，如果被积函数用积分下限处的函数值代替被积函数作积分，则有

()()()()()()11,,k k

t k k t k k k y t y t f t y t dt

y t hf t y t ++=+??≈+??

?

进而得到下式给出的递推算法—Euler 方法

()

()100

,1,2,,k k k k y y hf t y k n

y t y +=+==

2、Euler 方法的改进

改进的Euler 方法Euler 方法的精度高，其原因在于： 在推导Euler 方法时，我们是用待求解函数在一点处的变化率

()()

,k k f t y t （1）

代替在区间上的平均变化率；

而在推导改进的Euler 方法时，我们是用待求解函数在两点处变化率的平均值

()()()()111

,,2k k k k f t y t f t y t ++??

+?

?（2）

代替在区间上的平均变化率；

显然，通常式（2）比式（1）更接近于在区间上的平均变化率。由此，人们受到启发：适当地选取区间上函数多个点处的变化率，用它们加权平均值代替在区间上的平均变化率，近似解的精度应更高。

二、 Runge —Kutta 法：①二阶的Runge —Kutta 法②三阶的Runge —Kutta 法③四阶的

Runge —Kutta 法