同济大学数值分析matlab编程Word文档下载推荐.docx

同济大学数值分析matlab编程Word文档下载推荐.docx

《同济大学数值分析matlab编程Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《同济大学数值分析matlab编程Word文档下载推荐.docx（18页珍藏版）》请在冰豆网上搜索。

end

editf.m

functionv=f（x）

v=x.*log（x）-1;

editg.m

functionz=g（y）

z=y.^5+y-1;

[x1iter1]=gexianfa（'

f'

1,3,1e-10）

x1=

1.7632

iter1=

6

[x2iter2]=gexianfa（'

g'

0,1,1e-10）

x2=

0.7549

iter2=

8

3.使用GS迭代求解下述线性代数方程组：

editgsdiedai.m

function[xiter]=gsdiedai（A,x0,b,tol）

D=diag（diag（A））;

L=D-tril（A）;

U=D-triu（A）;

x=x0;

while（（norm（b-A*x）./norm（b））>

x0=x;

x=（D-L）\（U*x0+b）;

A=[521;

-142;

1-310];

b=[-12103]'

tol=1e-4;

x0=[000]'

[xiter]=gsdiedai（A,x0,b,tol）;

x

x=

-3.0910

1.2372

0.9802

iter

iter=

4.用四阶Range-kutta方法求解下述常微分方程初值问题（取步长h=0.01）

editksf2.m

functionv=ksf2（x,y）

v=y+exp（x）+x.*y;

a=1;

b=2;

h=0.01;

n=（b-a）./h;

x=[1:

0.01:

2];

y

（1）=2;

fori=2:

（n+1）

k1=h*ksf2（x（i-1）,y（i-1））;

k2=h*ksf2（x（i-1）+0.5*h,y（i-1）+0.5*k1）;

k3=h*ksf2（x（i-1）+0.5*h,y（i-1）+0.5*k2）;

k4=h*ksf2（x（i-1）+h,y（i-1）+k3）;

y（i）=y（i-1）+（k1+2*k2+2*k3+k4）./6;

y

调用函数方法

editRangekutta.m

function[xy]=Rangekutta（f,a,b,h,y0）

x=[a:

h:

b];

n=（b-a）/h;

y

（1）=y0;

k1=h*（feval（f,x（i-1）,y（i-1）））;

k2=h*（feval（f,x（i-1）+0.5*h,y（i-1）+0.5*k1））;

k3=h*（feval（f,x（i-1）+0.5*h,y（i-1）+0.5*k2））;

k4=h*（feval（f,x（i-1）+h,y（i-1）+k3））;

[xy]=Rangekutta（'

ksf2'

1,2,0.01,2）;

5.取

，请编写Matlab程序，分别用欧拉方法、改进欧拉方法在

上求解初值问题。

editEuler.m

function[xy]=Euler（f,a,b,h,y0）

n=（b-a）./h;

y（i）=y（i-1）+h*feval（f,x（i-1）,y（i-1））;

editgaijinEuler.m

function[xy]=gaijinEuler（f,a,b,h,y0）

y1=y（i-1）+h*feval（f,x（i-1）,y（i-1））;

y2=y（i-1）+h*feval（f,x（i）,y1）;

y（i）=（y1+y2）./2;

editksf3.m

functionv=ksf3（x,y）

v=x.^3-y./x;

[xy]=Euler（'

ksf3'

1,2,0.2,0.4）

1.00001.20001.40001.60001.80002.0000

y=

0.40000.52000.77891.21651.88362.8407

[xy]=gaijinEuler（'

0.40000.58950.92781.46152.24643.3466

6.请编写复合梯形积分公式的Matlab程序，计算下面积分的近似值，区间等分

。

编写辛普森积分公式的Matlab程序，计算下面积分的近似值，区间等分

、

edittixingjifen.m

functions=tixingjifen（f,a,b,n）

x=linspace（a,b,（n+1））;

y=zeros（1,length（x））;

y=feval（f,x）

h=（b-a）./n;

s=0.5*h*（y

（1）+2*sum（y（2:

n））+y（n+1））;

editsimpson.m

functionI=simpson（f,a,b,n）

h=（b-a）/n;

x=linspace（a,b,2*n+1）;

y=feval（f,x）;

I=（h/6）*（y

（1）+2*sum（y（3:

2:

2*n-1））+4*sum（y（2:

2*n））+y（2*n+1））;

editksf4.m

functionv=ksf4（x）

v=1./（x.^2+1）;

tixingjifen（'

ksf4'

0,1,20）

ans=

0.7853

simpson（'

0,1,10）

0.7854

editksf5.m

functionv=ksf5（x）

if（x==0）

v=1;

else

v=sin（x）./x;

（第二个函数‘ksf5’调用求积函数时，总显示有错误：

“NaN”，还没调试好。

见谅！

7.用

迭代方法对下面方程组求解，取初始向量

editJacobi.m

function[xiter]=Jacobi（A,x0,b,tol）

while（norm（A*x-b）/norm（b）>

x=D\（（L+U）*x0+b）;

A=[24-4;

333;

442];

b=[2-3-2]'

x0=[32-1]'

[x,iter]=Jacobi（A,x0,b,1e-4）

1

-1

3

8.用牛顿法求解方程

在

附近的根。

editNewton.m

function[xiter]=Newton（f,g,x0,tol）

x=x0-feval（f,x0）./feval（g,x0）;

editksf6.m

functionv=ksf6（x）

v=x*cos（x）+2;

editksg6.m

functionz=ksg（y）

[xiter]=Newton（'

ksf6'

'

ksg6'

2,1e-4）

2.4988

9.分别用改进乘幂法、反幂法计算矩阵A的按模最大特征值及其对应的特征向量、按模最小特征值及其对应的特征向量。

解

editep.m

function[t,x]=ep（A,x0,tol）

[tv0ti0]=max（abs（x0））;

lam0=x0（ti0）;

x0=x0./lam0;

x1=A*x0;

[tv1ti1]=max（abs（x1））;

lam1=x1（ti1）;

x1=x1./lam1;

while（abs（lam0-lam1）>

lam0=lam1;

[tv1ti1]=max（abs（x1））;

lam1=x1（ti1）;

t=lam1;

editfanep.m

function[t,x]=fanep（A,x0,tol）

x1=A\x0;

while（abs（1/lam0-1/lam1）>

t=1/lam1;

A=[126-6;

6162;

-6216];

x0=[10.5-0.5]'

[t,x]=ep（A,x0,tol）

t=

21.5440

1.0000

0.7953

-0.7953

[t,x]=fanep（A,x0,tol）

4.4560

-0.6287

0.6287

10.将积分区间n等分，用复合梯形求积公式计算定积分

比较不同

值时的误差（画出平面上log（n）-log（Error）图）．

editksf7.m

functionv=ksf7（x）

v=sqrt（1+x.^2）;

I=quad（'

ksf7'

1,3）;

n=1:

100;

fori=1:

100

x（i）=tixingjifen（'

1,3,i）;

error（i）=abs（I-x（i））;

plot（log10（n）,log10（error（n）））

11.用

editsor.m

function[x,iter]=sor（A,x0,b,omega,tol）

x=（D-omega*L）\（omega*b+（1-omega）*D*x+omega*U*x）;

A=[2-10;

-13-1;

0-12];

b=[18-5]'

x0=[10-1]'

[x,iter]=sor（A,x0,b,1.1,1e-4）

1.9999

3.0000

-1.0000

5

其他

A=[1234;

5678];

A

A=

1234

5678

[mn]=size（A）

m=

2

n=

4

x=[12345];

12345

length（x）

A=[234;

119;

12-6];

234

119

12-6

[LU]=lu（A）;

L

L=

1.000000

0.50001.00000

0.5000-1.00001.0000

U

U=

2.00003.00004.0000

0-0.50007.0000

00-1.0000

x=linspace（0,1,11）;

x'

0

0.1000

0.2000

0.3000

0.4000

0.5000

0.6000

0.7000

0.8000

0.9000

1.0000

x=[1234];

y=[6111827];

p=polyfit（x,y,2）

p=

1.00002.00003.0000

diag（ones（4,1）,1）

01000

00100

00010

00001

00000

1牛顿法解方程组根

Functions=NewtonIterate（x,eps）

%Newton迭代法求解非线性方程组的解

%x为迭代初值，eps为允许误差

Ifnargin==1

eps=1.0e-6;

elseifnargin<

1

return

x1=fx1（x）;

%非线性方程组函数fx1.m

x2=-dfx1（x）;

%非线性方程组的导数函数dfx1.m

x0=x2\x1'

whilenorm（x0）>

=eps

x=x0'

+x;

s=x0'

functiony=fx1（x）%函数1

y=x^3-2*x-5

functiony=dfx1（x）%函数2

y=3*x^2-2

%命令行

x=NewtonIterate（2.5,eps）

x=1307/624

2拉格朗日插值法

functionyh=lagrange（x,y,xh）

n=length（x）;

m=length（xh）;

yh=zeros（1,m）;

c1=ones（n-1,1）;

c2=ones（1,m）;

n

xp=x（[1:

i-1i+1:

n]）;

yh=yh+y（i）*prod（（c1*xh-xp'

*c2）./（x（i）-xp'

*c2））;

3牛顿插值法

function[p,q]=chashang（x,y）

p（:

1）=x;

2）=y;

forj=3:

n+1

p（1:

n+2-j,j）=diff（p（1:

n+3-j,j-1））./（x（j-1:

n）-x（1:

n+2-j））;

q=p（1,2:

n+1）'

4改进乘幂法

function[t,y]=eigIPower（a,xinit,ep）

v0=xinit;

[t,ti]=max（abs（v0））;

lam0=v0（ti）;

u0=v0/lam0;

flag=0;

while（flag==0）

v1=a*u0;

[tv,ti]=max（abs（v1））;

lam1=v1（ti）;

u0=v1/lam1;

err=abs（lam0-lam1）;

if（err<

=eo）

flag=1;

y=u0;

5高斯赛德尔迭代法求解方程组

function[x,iter]=gs（A,b,tol）

x=zeros（size（b））;

foriter=1:

500

x=（D-L）\（b+U*x）;

error=norm（b-A*x）/norm（b）;

if（error<

break;

6复合梯形公式

7复合辛普森公式

8欧拉公式求解微分方程

function[x,y]=odeEuler（f,y0,a,b,n）

x=a:

b;

n,

y（i+1）=y（i）+h*feval（f,x（i）,y（i））;