matlab程序设计实践

1.编程实现以下科学计算算法，并举一例应用之。（参考书籍《精通

MATLAB科学计算》，王正林等著，电子工业出版社，2009年）“牛顿下山法求非线性方程组解”

解：

算法说明：

牛顿下山法的迭代公式为：

X n+1= x n -ω(F(x n))-1F(x n)

ω的取值范围为0<ω≤1,为了保证收敛，还要要求ω的取值使得：

‖F(x n+1)‖<‖F(x n)‖

可以用逐次减半法来确定ω。为了减少计算量，还可以用差商来代替偏导数。在MATLAB中编程实现的非线性方程组的牛顿下山法的函数为：mulDNewton。

功能：用牛顿下山法求非线性方程组的一个解。

调用格式：[r,n]=mulDNewton(x0,eps)。

其中，x0为初始迭代向量；

eps为迭代精度；

r为求出的解向量；

n为迭代步数。

牛顿下山法的MATLAB代码如下：

function [r,n]=mulDNewton(x0,eps)

%牛顿下山法求非线性方程组的一个解

%初始迭代向量：x0

%迭代精度：eps

%解向量：r

%迭代步数：n

if nargin==1

eps=1.0e-4;

end

r=x0-myf(x0)/dmyf(x0);

n=1;

tol=1;

while tol>eps

x0=r;

ttol=1;

w=1;

F1=norm(myf(x0));

while ttol>=0 %下面的循环是选取下山因子w的过程

r=x0-w*myf(x0)/dmyf(x0); %核心的迭代公式

ttol=norm(myf(r))-F1;

w=w/2;

end

tol=norm(r-x0);

n=n+1;

if (n>100000) %迭代步数控制

disp('迭代步数太多，可能不收敛！');

return;

end

end

实例：牛顿下山法求解非线性方程组应用实例。采用牛顿下山法求下面方程组的一个根。

其初始迭代值取（0,0）。

解：首先建立myf.m函数文件，输入以下内容：

function f=myf(x)

f(1)=0.5*sin(x(1))+0.1*cos(x(2)*x(1))-x(1);

f(2)=0.5*cos(x(1))-0.1*sin(x(2))-x(2);

f=[f(1) f(2)];

再建立dmyf.m导数的雅克比矩阵，输入以下内容：

function df=dmyf(x)

df=[0.5*cos(x(1))-0.1*x(2)*sin(x(2)*x(1))-1 -0.1*x(1)*sin(x(2)*x(1));

-0.5*sin(x(1)) -0.1*cos(x(2))-1];

然后，在MATLAB命令窗口中输入：

>> [r,n]=mulDNewton([0,0])

输出计算结果为：

r =

0.1979 0.4470

n =

5

由计算结果可知，初始迭代值取（0,0）时，用5步迭代得到了方程组的一组解（0.1979,0.4470），牛顿下山法也是一种快速且有效的求解非线性方程组的方法。

流程图：

⑴源程序

⑵例题

2．编程解决以下科学计算和工程实际问题

⑴

2-4在平炉炼钢过程中，由于矿石及炉气的氧化作用，铁水中的总含碳量在不断降低，一炉钢在冶炼初期（熔化期）中总的去碳量y与所加的两种矿石（天然矿石与烧结矿石）的量x1,x2及熔化时间x3有关（熔化时间愈长则去碳量愈多）。经实测，某号平炉的49炉数据如表2.6所示，试求出y与x1,x2,x3之间的线性回归方程。

解：

在MATLAB命令窗口中输入：

>> x1 =[2 7 5 12 1 3 2 6 7 0 3 0 8 6 0 3 7 16 6 0 9 4 0 9 2 9 12 6 12 0 5 4 0 6 4 10 4 5 9 6 5 5 8 2 7 4 10 3 4];

>> x2 =[18 9 14 3 20 12 17 5 8 23 16 18 4 14 21 14 12 0 16 15 0 6 17 0

16 6 5 13 7 24 12 15 20 16 17 4 14 13 8 13 8 11 6 13 8 10 5 17 15];

>> x3 =[50 40 46 43 64 40 64 39 37 55 60 49 50 51 51 51 56 48 46 52 40 32 47 44 39 39 61 41 47 61 37 48 45 42 48 48 36 36 51 54 100 44 63 55 50 45 40 64 72];

>> y=[4.3302 3.6485 4.4830 5.5468 5.4970 3.1125 5.1182 3.8759 4.6700

4.9536

5.0060 5.2701 5.8772 5.4849 4.5960 5.6645

6.0795 3.2194

5.8076 4.7306 4.6805 3.1272 2.6104 3.7174 3.8946 2.7066 5.6314 5.8152 5.1302 5.3910 4.4533 4.6569 4.5212 4.8650 5.3566 4.6096 2.3815 3.8746 4.5919 5.1588 5.4373 3.9960 4.3970 4.0622 2.2905 4.7115 4.5310 5.3637

6.0771];

>> y=y';

>> x=[x1',x2',x3',ones(49,1)];

>> regress(y,x)

ans =

0.1504

0.1014

0.0370

0.7753

得到：

0.1504y=0.1014x1+0.0370x2+0.7753x3

流程图：