matlab数值微积分与方程数值求解

电子一班王申江

实验九数值微积分与方程数值求解

一、实验目的

1、掌握求数值导数和数值积分的方法

2、掌握代数方程数值求解的方法

3、掌握常微分方程数值求解的方法

二、实验内容

1、求函数在指定点的数值导数。

()

23

2

123,1,2,3 026

x x x

f x x x x

x

==

>>syms x

>>f=[x x^2 x^3;1 2*x 3*x^2;0 2 6*x]; >>F=det(f)

F=2*x^3

>>h=0.1

>>x=[0:h:4];

>>f=2*x^3;

>>[dy,dx]=diff_ctr(f,h,1);

>>y1=dy(dx==1)

y1=6.0000

>>y2=dy(dx==2)

y2=24.0000

>>y3=dy(dx==3) y3=54.0000

2、用数值方法求定积分。

（1）

210I π

=⎰的近似值 a=inline('sqrt(cos(t.^2)+4*sin((2*t).^2)+1)');

I=quadl(a,0,2*pi)

I =

6.7992 + 3.1526i

（2）()1

202ln 11x I dx x

+=+⎰ b=inline('log(1+x)./(1+x.^2)');

I=quadl(b,0,1)

I =

0.2722

3、分别用3种不同的数值方法解线性方程组。

6525494133422139211

x y z u x y z u x y z u x y u +-+=-⎧⎪-+-=⎪⎨++-=⎪⎪-+=⎩

A=[6,5,-2,5;9,-1,4,-1;3,4,2,-2;3,-9,0,2];

b=[-4,13,1,11]';

x=A\b

x =

0.6667

-1.0000

1.5000

-0.0000

4、求非齐次线性方程组的通解。

123412341

234273635224

9472x x x x x x x x x x x x +++=⎧⎪+++=⎨⎪+++=⎩

A=[2 7 3 1;3 5 2 2 ;9 4 1 7]

A =

2 7

3 1

3 5 2 2

9 4 1 7

>> b=[6 4 2]'

b =

6

4

2

>> rank(A)

ans = 2

>> rank([A,b])

ans =

2

>> x=pinv(A)*b

x =

-0.0571

0.7429

0.3429 -0.1143

5、求代数方程的数值解。

（1）03sin 0 1.5x x x e x +-==在附近的根

function f=fz(x)

f=3*x-sin(x*pi/180)-exp(x)

fzero('fz',-1.5)

ans =

0.6287

（2）在给定的初值0001,1,1x y z ===下，求方程组的数值解。

2sin ln 703231050y x y z x z x y z ⎧++-=⎪+-+=⎨⎪++-=⎩

function F=myfun (X) x=X(1);

y=X(2);

z=X(3);

F(1)=sin(x)+y^2+log(z)-7;

F(2)=3*x+2^y-3*z+1;

F(3)=x+y+z-5;

X=fsolve('myfun',[1,1,1],optimset('Display', 'off')) X = 0.2136 2.4359 2.3505

6、求函数在指定区间的极值。

（1）()3cos log x

x x x x f x e ++=在（0，1）内的最小值。 function f=f(x)

f=(x^3+cos(x)+x+x*log(x))/exp(x);

[x,fval]=fminbnd(f,0,1)

（2）()321,21121222410f x x x x x x x x =+-+在[0，

0]附近的最小值点和最小值。

function fxy(u)

x=u(1);

y=u(2);

f=2.*x.^2-6.*x.*y+y.^2;

[U,fmin]=fminsearch('fxy',[0,0])

7、求微分方程的数值解，并绘制解的曲线。

()()22'500000xd y dy y dx dx y y ⎧-+=⎪⎪⎪=⎨⎪=⎪⎪⎩

function dy=vdp5(x,y)

dy=zeros(2,1)

dy(1)=y(2)

dy(2)=5*y(2)/x-y(1)/x

[X,Y]=ode23('vdp5',[0,100],[0,0]);

plot(X,Y(:,1),'-')

8、求微分方程组的数值解，并绘制解的曲线。()()()'1'213'3121230.5100,01,01y yy y y y y y y y y y ⎧=⎪=-⎪⎨=-⎪⎪===⎩

function dy=rigid(t,y)

dy=zeros(3,1);

dy(1)=y(2)*y(3);

dy(3)=-0.51*y(1)*y(2);

[T,Y]=ode45('rigid',[0,20],[0,1,1]);

plot(T,Y(:,1),'-',T,Y(:,2),'+',T,Y(:,3),'*')