数值计算常见编程汇总20教材.docx

数值计算常见编程汇总20教材.docx

《数值计算常见编程汇总20教材.docx》由会员分享，可在线阅读，更多相关《数值计算常见编程汇总20教材.docx（26页珍藏版）》请在冰豆网上搜索。

数值计算常见编程汇总20教材

1.秦九韶法

利用秦九韶算法求多项式p（x）=x^5-3x^4+4x^2-x+1在x=3时的值。

x=input（'请输入x：

'）;

n=input（'请输入多项式最高幂数n:

'）;

a=input（'请按照升幂顺序写出系数矩阵：

'）;

v0=a（n+1）;

fork=1:

n

v1=v0*x+a（n+1-k）;

v0=v1;

end

v=v0

请输入x：

3

请输入多项式最高幂数n:

5

请按照升幂顺序写出系数矩阵：

[1-140-31]

v=34

2.二分法

用二分法求方程x^3-x-1=0在[1,2]内的近似根，要求误差不超过10-3

a=1;

b=2;

fa=a^3-a-1;

fb=b^3-b-1;

c=（a+b）/2;

fc=c^3-c-1;

iffa*fb>0,break,end

whileabs（fc）>0.001

c=（a+b）/2;

fc=c^3-c-1;

iffa*fc<0;

b=c;

fb=fc;

else

a=c;

fa=fc;

end

end

formatlong

fx=fc,x=c

fx=-4.659488331526518e-05

x=1.324707031250000

3.拉格朗日差值法：

已知点（10）、（27）、（316）求x=4的值。

x=input（'以行向量形式输入x的取值'）;

y=input（'以行向量形式输入y的取值'）;

n=input（'多项式的次数'）;

x1=input（'计算在x处的函数值,x1='）;

fork=1:

（n+1）

l（k）=1;

end

fork=1:

（n+1）

fori=1:

k-1

l（k）=l（k）*（x1-x（i））/（x（k）-x（i））;

end

fori=（k+1）:

（n+1）

l（k）=l（k）*（x1-x（i））/（x（k）-x（i））;

end

end

p=0;

fork=1:

n+1

p=p+l（k）*y（k）;

end

disp（p）;

以行向量形式输入x的取值[123]

以行向量形式输入y的取值[0716]

多项式的次数2

计算在x处的函数值,x1=4

27

4.分段插值法：

已知点（00.9）、（43.4）、（65.8）（88.9）（1012）求x=7.5的值。

x=[046810];

y=[0.93.45.88.912];

n=length（x）;

x0=input（'输入待求x0的值：

'）;

fori=1:

n-1

ifx0

（1）||x0>x（n）

disp（'不符合条件，无法进行分段插值'）;

break;

end

ifx（i）<=x0&&x0

f=y（i）+（y（i+1）-y（i））/（x（i+1）-x（i））*（x0-x（i））;

end

end

F

输入待求x0的值：

7.5

f=8.1250

5.牛顿差值

已知点（23），（37），（59），（812）求x=7的值

x0=2;x1=3;x2=5;x3=8;

y0=3;y1=7;y2=9;y3=12;

x=7;

f11=（y1-y0）/（x1-x0）;

f12=（y2-y0）/（x2-x0）;

f13=（y3-y0）/（x3-x0）;

f22=（f12-f11）/（x2-x1）;

f23=（f13-f11）/（x3-x1）;

f33=（f23-f22）/（x3-x2）;

N=y0+f11*（x-x0）+f22*（x-x1）*（x-x0）+f33*（x-x2）*（x-x1）*（x-x0）;

p=N

p=

9.6667

6.分段三次插值：

已知：

p（0）=0,p

（1）=3,p

（2）=12,p（3）=33,p，（0）=2，p，

（1）=5，p，

（2）=14，p，（3）=29，用分段三次插值法求解p（1.2）

clc;

x=input（'已知x的取值'）;

y=input（'已知y的取值'）;

yy=input（'已知y的一阶导数的取值'）;

n=input（‘请输入x值的个数：

’）;

x0=input（'x0的取值'）;y0=0;

fori=1:

n-1

if（x0

（1）||x0>x（n））

disp（'输入的值不在所给范围内'）;break;

end

if（x（i）<=x0&&x0

h=x（i+1）-x（i）;

t=（x0-x（i））/h;y0=（t-1）^2*（2*t+1）*y（i）+t^2*（3-2*t）*y（i+1）+h*t*（t-1）^2*yy（i）+h*t^2*（t-1）*yy（i+1）;

end

end

disp（y0）;

已知x的取值[0,1,2,3]

已知y的取值[0,3,12,33]

已知y的一阶导数的取值[2,5,14,29]

请输入x值的个数：

4

x0的取值1.2

4.1280

7.牛顿插值公式

已知：

p（0）=0,p

（1）=3,p

（2）=12,p（3）=33,用牛顿插值法求解p（1.2）

clc;

t=input（'以行向量形式输入x的取值'）;

y=input（'以行向量形式输入y的取值'）;

n=input（'多项式的次数'）;

x=input（'计算在x处的函数值,x='）;

fork=2:

（n+1）

f（k）=0;

end

fork=2:

（n+1）

fori=1:

k

u=1;

forj=1:

（i-1）

u=u*（t（i）-t（j））;

end

forj=（i+1）:

k

u=u*（t（i）-t（j））;

end

f（k）=f（k）+y（i）/u;，，

end

end

p=y

（1）;

fork=2:

（n+1）

v=1;

fori=1:

（k-1）

v=v*（x-t（i））;

end

p=p+f（k）*v;

end

disp（p）;

以行向量形式输入x的取值[0,1,2,3];

以行向量形式输入y的取值[0,3,12,33]

多项式的次数3

计算在x处的函数值,x=1.2

4.1280

8.复化梯形积分公式

已知f（x）=3x^2+2*x,将区间[0,3]划分为3段，用复化梯形积分公式求函数f（x）在区间[0,3]积分。

方法1

clc;

f=input（'请输入函数形式','s'）;

f=inline（f）;

a=input（'积分下限a：

'）;

b=input（'积分上限b：

'）;

n=input（'需要划分的段数：

'）;

h=（b-a）/n;

x

（1）=a;

x（n+1）=b;

s=f（x

（1））+f（x（n+1））;

fori=2:

n

x（i）=x（i-1）+h;

s=s+2*f（x（i））;

end

s=s*h/2;

disp（s）;

请输入函数形式3*x^2+2*x

积分下限a：

0

积分上限b：

3

需要划分的段数：

3

37.5000

方法2

clc;

a=input（'积分下限a：

'）;

b=input（'积分上限b：

'）;

n=input（'需要划分的段数：

'）;

h=（b-a）/n;

x

（1）=a;

x（n+1）=b;

fori=1:

n+1

x（i）=x

（1）+（i-1）*h;

f（i）=3*x（i）^2+2*x（i）;

end

s=f

（1）+f（n+1）;

fori=2:

n

s=s+2*f（i）;

end

s=s*h/2;

disp（s）;

积分下限a：

0

积分上限b：

3

需要划分的段数：

3

37.5000

9.复化辛普生积分公式

已知f（x）=3x^2+2*x,将区间[0,3]划分为3段，用复化辛普生积分公式求函数f（x）在区间[0,3]积分。

方法1

clc;

f=input（'请输入函数形式','s'）;

f=inline（f）;

a=input（'积分下限a：

'）;

b=input（'积分上限b：

'）;

n=input（'需要划分的段数：

'）;

h=（b-a）/n;

x

（1）=a;

x（n+1）=b;

s=f（x

（1））+f（x（n+1））;

fori=2:

n

x（i）=x（i-1）+h;

s=s+4*f（（x（i-1）+x（i））/2）+2*f（x（i））;

end

s=s+4*f（（x（n）+x（n+1））/2）;

s=s*h/6;

disp（s）;

请输入函数形式3*x^2+2*x

积分下限a：

0

积分上限b：

3

需要划分的段数：

3

36

方法2

clc;

a=input（'积分下限a：

'）;

b=input（'积分上限b：

'）;

n=input（'需要划分的段数：

'）;

h=（b-a）/n;

x

（1）=a;

x（n+1）=b;

fori=1:

n+1

x（i）=x

（1）+（i-1）*h;

f（i）=3*x（i）^2+2*x（i）;

end

fori=1:

n

y（i）=（x（i）+x（i+1））/2

l（i）=3*y（i）^2+2*y（i）;

end

s=f（n+1）-f

（1）;

fori=1:

n

s=s+4*l（i）+2*f（i）;

end

s=s*h/6;

disp（s）;

积分下限a：

0

积分上限b：

3

需要划分的段数：

3

36

10.两点高斯公式

已知f（x）=3x^2+2*x,将区间[0,3]划分为3段，用两点高斯积分公式求函数f（x）在区间[0,3]积分。

clc;

f=input（'请输入函数形式','s'）;

f=inline（f）;

a=input（'积分下限a：

'）;

b=input（'积分上限b：

'）;

n=input（'需要划分的段数n：

'）;

h=（b-a）/n;

t=2*3^0.5;

x

（1）=a;

s=0;

fori=2:

n+1

x（i）=x（i-1）+h;

c=（x（i）+x（i-1））/2;

s=s+f（c-h/t）+f（c+h/t）;

end

s=s*h/2;

disp（s）;

请输入函数形式3*x^2+2*x

积分下限a：

0

积分上限b：

3

需要划分的段数n：

3

36

11.复化高斯公式

已知f（x）=3x^2+2*x,将区间[0,3]划分为3段，用两点高斯积分公式求函数f（x）在区间[0,3]积分.

程序

clc;

f=input（'请输入函数形式','s'）;

f=inline（f）;

a=input（'积分下限a：

'）;

b=input（'积分上限b：

'）;

n=input（'需要划分的段数n：

'）;

h=（b-a）/n;

x

（1）=a;

x（n+1）=b;

s=0;

fori=2:

n

x（i）=x（i-1）+h;

end

fork=1:

n

s=s+f（x（k）+h/2+h/（2*3^0.5））+f（x（k）+h/2-h/（2*3^0.5））;

end

s=s*h/2;

disp（s）;

程序运行结果：

请输入函数形式3*x^2+2*x

积分下限a：

0

积分上限b：

3

需要划分的段数n：

3

36

12.导数公式

已知：

p（0）=0,p

（1）=3,p

（2）=12,分别用向前求导、中点求导、向后求导求p,（0）、p,

（1）、p,

（2）

clc;

x0=input（'请输入x0的值：

'）;

y0=input（'请输入y0的值：

'）;

x1=input（'请输入x1的值：

'）;

y1=input（'请输入y1的值：

'）;

x2=input（'请输入x2的值：

'）;

y2=input（'请输入y2的值：

'）;

h=input（'请输入h步长的值：

'）;

f=（y1–y0）/h;

g=（y2-y0）/（2*h）;

k=（y2-y1）/h;

disp（f）;

disp（g）;

disp（k）;

请输入x0的值：

0

请输入y0的值：

0

请输入x1的值：

1

请输入y1的值：

3

请输入x2的值：

2

请输入y2的值：

12

请输入h步长的值：

1

3

6

9

13.一步欧拉方法

已知：

y’=x+xy，y（0）=1,h=0.2,用一步欧拉方法求解y（0.8）

方法1、

clearall;

clc;

f=input（'请输入导数的函数形式','s'）;

f=inline（f）;

x0=input（'x0的初始值'）;

y0=input（'y0的初始值'）;

h=input（'步长'）;

x=input（'x的值'）;

while（x0

x1=x0+h;

y1=y0+h*f（x0,y0）;

x0=x1;

y0=y1;

end

disp（y0）;

请输入导数的函数形式x+x*y

x0的初始值0

y0的初始值1

步长0.2

x的值0.8

1.5160

方法2、

clearall;

clc;

x1=input（'x1的初始值:

'）;

y1=input（'y1的初始值:

'）;

h=input（'步长:

'）;

x=input（'x的值:

'）;

n=（x-x1）/h;

fori=1:

n+1

y（i）=1

end

fori=1:

n

x（i）=x1+（i-1）*h;

y（i+1）=y（i）+h*（x（i）+x（i）*y（i））;

end

disp（y（5））

x1的初始值:

0

y1的初始值:

1

步长:

0.2

x的值:

0.8

1.5160

14.二步欧拉方法

已知：

y’=x+xy，y（0）=1,h=0.2,用二步欧拉方法求解y（0.8）

方法1

clearall;

clc;

f=input（'请输入导数的函数形式','s'）;

f=inline（f）;

x0=input（'x0的初始值：

'）;

y0=input（'y0的初始值：

'）;

h=input（'步长h：

'）;

x=input（'x的值：

'）;

x1=x0+h;

y1=y0+h*f（x0,y0）;

while（x1

x2=x1+h;

y2=y0+2*h*f（x1,y1）;

y0=y1;y1=y2;

x0=x1;x1=x2;

end

disp（y1）;

请输入导数的函数形式x+x*y

x0的初始值：

0

y0的初始值：

1

步长h：

0.2

x的值：

0.8

1.7229

方法2：

clearall;

clc;

x0=input（'x0的初始值：

'）;

y0=input（'y0的初始值：

'）;

h=input（'步长h：

'）;

x=input（'x的值：

'）;

n=（x-x0）/h;

fori=1:

n+1

y（i）=1;

end

y

（2）=y0+h*（x0+x0*y0）;

fori=1:

n-1

x（i+1）=x0+i*h;

y（i+2）=y（i）+2*h*（x（i+1）+x（i+1）*y（i+1））;

end

disp（y（5））

x0的初始值：

0

y0的初始值：

1

步长h：

0.2

x的值：

0.8

1.7229

15.隐式欧拉方法及其改进方法

已知：

y’=x+xy，y（0）=1,h=0.2,用二步欧拉方法求解y（0.8）

clearall;

clc;

f=input（'请输入导数的函数形式','s'）;

f=inline（f）;

x0=input（'x0的初始值：

'）;

y0=input（'y0的初始值：

'）;

h=input（'步长h：

'）;

xn=input（'终值：

'）;

while（x0

x1=x0+h;

y1=y0+h*f（x0,y0）;

y0=y0+h/2*（f（x0,y0）+f（x1,y1））;

x0=x1;

end

disp（y1）;

请输入导数的函数形式x+x*y

x0的初始值：

0

y0的初始值：

1

步长h：

0.2

终值：

0.8

1.6797

16.方程组差分法

初值y（0）=1,z（0）=2,取步长h=0.2,用方程组差方法求得区间[0,1]之间的各值。

方法1：

clearall

clc

f=input（'输入y的导数函数：

'）;

g=input（'输入z的导数函数：

'）;

h=input（'输入步长取值：

'）;

a=input（'输入差分区间的下边界：

'）;

b=input（'输入差分区间的上边界：

'）;

N=（b-a）/h;

y=input（'输入y0的取值：

'）;

z=input（'输入z0的取值'）;

fork=1:

N

x=a+（k-1）*h;

c=y;

y=y+h*f（x,y,z）;Y（k）=y;

z=z+h*g（x,y,z）;Z（k）=z;

end

formatshort

Y

Z

输入y的导数函数：

inline（'3*x^2+z','x','y','z'）

输入z的导数函数：

inline（'y+z','x','y','z'）

输入步长取值：

0.2

输入差分区间的下边界：

0

输入差分区间的上边界：

1

输入y的取值：

1

输入z的取值：

2

Y=1.40001.94402.72003.82985.3951

Z=2.60003.40004.46885.90667.8538

方法2：

h=input（'输入步长取值：

'）;

a=input（'输入差分区间的下边界：

'）;

b=input（'输入差分区间的上边界：

'）;

N=（b-a）/h;

y

（1）=input（'输入y0的取值：

'）;

z

（1）=input（'输入z0的取值'）;

fork=1:

N

x（k）=a+（k-1）*h;

y（k+1）=y（k）+h*（3*x（k）^2+z（k））;

z（k+1）=z（k）+h*（y（k）+z（k））;

end

formatshort

disp（y）

disp（z）

输入步长取值：

0.2

输入差分区间的下边界：

0

输入差分区间的上边界：

1

输入y0的取值：

1.0

输入z0的取值:

2.0

1.00001.40001.94402.72003.82985.3951

2.00002.60003.40004.46885.90667.8538

17.一般迭代法：

求方程exp（-x）=Ax的全部根；初值自己确定，A=5

程序：

方法1

clearall;

clc;

f=inline（'exp（-x）/5','x'）;

esp=10^（-4）;

x0=0;

x1=f（x0）;

while（abs（x1-x0）>esp）

x2=f（x1）;

x0=x1;

x1=x2;

end

formatshorte

x1

程序运行结果：

x1=1.6892e-01

方法2

clearall;

clc;

n=input（‘请输入迭代次数：

’）

x0=0;

fori=1:

n

x1=0.2*exp（-x0）;

x0=x1;

end

formatshorte

x0

请输入迭代次数：

5

x0=

1.6894e-01

18.雅克比迭代法

用雅克比迭代法求解方程

a=input（'输入加工后方程组的系数矩阵a'）;

b=input（'输入加工后b的取值'）;

x=input（'输入方程组的解的迭代初值'）;

n=numel（b）;

N=input（'输入迭代的次数'）;

fori=1:

N

forj=1:

n

t（j）=b（j）;

fork=1:

n

t（j）=t（j）+a（j,k）*x（k）;

end

end

x=t;

end

disp（x）;

输入加工后方程组的系数矩阵a[00-0.10.5;-1/803/80;3/81/401/8;-1/72/7-2/70]

输入加工后b的取值[-0.711/8-23/817/7]

输入方程组的解的迭代初值[0000]

输入迭代的次数15

1.00000.5000-2.00003.0000

19.高斯--塞德尔迭代法

用高斯-----赛德尔迭代法求解方程

clearall;

clc;

a=input（'输入加工后方程组的系数矩阵a'）;

b=input（'输入加工后b的取值'）;

x=input（'输入方程组的解的迭代初值'）;

n=numel（b）;

N=input（'输入迭代的次数'）;

fori=1:

N

forj=1:

n

x（j）=b（j）;

fork=1:

n

x（j）=x（j）+a（j,k）*x（k）;

end

end

end

disp（x）;

输入加工后方程组的系数矩阵a[00-0.10.5;-1/803/80;3/81/401/8;-1/72/7-2/70]

输入加工后b的取值[-0.711/8-23/817/7]

输入方程组的解的迭代初值[0000]

输入迭代的次数10

1.00000.5000-2.00003.0000

20.追赶法

用追赶法求解方程

clearall;

clearall;

a=input（'输入方程组的系数a'）;

b=input（'输入方程组的系数b'）;

c=input（'输入方程组的系数c'）;

f=input（'输入方程组的系数f'）;

n=input（'输入方程组的未知数个数'）;

u

（1）=c

（1）/b

（1）;

y

（1）=f

（1）/b

（1）;

fork=2:

n

u（k）=c（k）/（b（k）-u（k-1）*a（k））;

y（k）=（f（k）-y（k-1）*a（k））/（b（k）-u（k-1）*a（k））;

end

x（n）=y（n）;

fori=1:

n-1

x（n-i）=y（n-i）-u（n-i）*x（n-i+1）;

end

disp（x）

输入方程组的系数a[0-1-1-3]

输入方程组的系数b[2325]

输入方程组的系数c[-1-2-40]

输入方程组的系数f[6101]

输入方程组的未知数个数4

30-2-1