matlab数学建模实例

第四周3.

中的三个根。

，在求8] [0,041.76938.7911.1-)(2

3=-+=x x x x f function y=mj()for x0=0:0.01:8

x1=x0^3-11.1*x0^2+38.79*x0-41.769;if (abs(x1)<1.0e-8)x0

end

end

4.分别用简单迭代法、埃特金法、牛顿法求解方程，并比较收敛性与收敛速度（ε分别取10-3、10-5、10-8）。

简单迭代法：

function y=jddd(x0)

x1=(20+10*x0-2*x0^2-x0^3)/20;k=1;

while (abs(x1-x0)>=1.0e-3)

x0=x1;

x1=(20+10*x0-2*x0^2-x0^3)/20;k=k+1;end x1k

埃特金法：

function y=etj(x0)

x1=(20-2*x0^2-x0^3)/10;x2=(20-2*x1^2-x1^3)/10;

x3=x2-(x2-x1)^2/(x2-2*x1+x0);k=1;

while (abs(x3-x0)>=1.0e-3)

x0=x3;

x1=(20-2*x0^2-x0^3)/10;x2=(20-2*x1^2-x1^3)/10;

x3=x2-(x2-x1)^2/(x2-2*x1+x0);k=k+1;end

2

,020102)(023==-++=x x x x x f

x3

k

牛顿法：

function y=newton(x0)

x1=x0-fc(x0)/df(x0);

k=1;

while(abs(x1-x0)>=1.0e-3)

x0=x1;

x1=x0-fc(x0)/df(x0);k=k+1;

end

x1

k

function y=fc(x)

y=x^3+2*x^2+10*x-20;

function y=df(x)

y=3*x^2+4*x+10;

第六周

1.解例6-4（p77)的方程组，分别采用消去法（矩阵分解）、Jacobi迭代法、Seidel迭代法、松弛法求解，并比较收敛速度。

消去法：

x=a\d

或

[L,U]=lu(a);

x=inv(U)inv(L)d

Jacobi迭代法：

function s=jacobi(a,d,x0)

D=diag(diag(a));

U=-triu(a,1);

L=-tril(a,-1);

C=inv(D);

B=C*(L+U);

G=C*d;

s=B*x0+G;

n=1;

while norm(s-x0)>=1.0e-8

x0=s;

s=B*x0+G;

n=n+1;

end

n

Seidel迭代法：

function s=seidel(a,d,x0)

D=diag(diag(a));

U=-triu(a,1);

L=-tril(a,-1);

C=inv(D-L);

B=C*U;

G=C*d;

s=B*x0+G;

n=1;

while norm(s-x0)>=1.0e-5

x0=s;

s=B*x0+G;

n=n+1;

end

n

松弛法：

function s=loose(a,d,x0,w)

D=diag(diag(a));

U=-triu(a,1);

L=-tril(a,-1);

C=inv(D-w*L);

B=C*((1-w)*D+w*U);

G=w*C*d;

s=B*x0+G;

n=1;

while norm(s-x0)>=1.0e-8

x0=s;

s=B*x0+G;

n=n+1;

end

n

2.练习MATLAB的常用矩阵语句，就龙格现象函数（p88）练习插值语句interp,spline，并比较。

3.测得血液中某药物浓度随时间的变化值为:

t(h)0.250.5 1.0 1.5 2.0 3.0 4.0 6.08.010.0

19.3018.1515.3614.1012.899.327.55 5.24 3.86 2.88 C(mg

/L)

求t=0.45, 1.75, 5.0, 6.0时的浓度C.

分别用n=4,5,9的拉格朗日插值计算；并用样条函数插值计算,并比较结果。

拉格朗日插值:

function s=lagr(n)

x=[0.250.51.01.52.03.04.06.08.010.0];

y=[19.3018.1515.3614.1012.899.327.555.243.862.88];

x0=[0.451.755.06.0];

m=length(x0);

for i=1:m

D=abs(x-x0(i));

I=1;

while I<=n+1

for a=1:length(x)

if D(a)==min(D)

c(I)=a;

D(a)=max(D)+1;

break

end

end

I=I+1;

end

b=sort(c);

z=x0(i);

t=0.0;

for k=1:length(b)

u=1.0;

for j=1:length(b)

if j~=k

u=u*(z-x(b(j)))/(x(b(k))-x(b(j)));

end

end

t=t+u*y(b(k));

end

s(i)=t;

end