Matlab求解非线性方程工程问题的作业2

题目二：超市有A、B、C三种果篮，同时装有苹果、香蕉、橘子三种水果。

现已知在A种果篮中三种水果质量比为2:2:1，B种果篮中为1:2:3，C种果篮中为3:1:1，超市制作一套A、
B、C三种果篮样品需消耗苹果6千克、香蕉6千克、橘子7.5千克。

求三种果篮各重几许？
一、列方程组
设A、B、C三种果篮分别重X、Y、Z千克
2X/5+Y/6+3Z/5=6
2X/5+2Y/6+Z/5=6
X/5+3Y/6+Z/5=7.5
二、Gauss法直接求解
1、输入高斯法程序，保存为Gauss.m文件
function [x,det,index]=Gauss(A,b)
[n,m]=size(A);nb=length(b);
if n~=m
error('The rows and columns of matrix A must be equal!');
return;
end
if m~=nb
error('The columns of A must be equal the length of b!');
return;
end
index=1;det=1;x=zeros(n,1);
for k=1:n-1
a_max=0;
for i=k:n
if abs(A(i,k))>a_max
a_max=abs(A(i,k));r=i;
end
end
if a_max<1e-10
index=0;return;
end
if r>k
for j=k:n
z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;
end
z=b(k);b(k)=b(r);b(r)=z;det=-det;
end
for i=k+1:n
m=A(i,k)/A(k,k);
for j=k+1:n
A(i,j)=A(i,j)-m*A(k,j);
end
b(i)=b(i)-m*b(k);
end
det=det*A(k,k);
end
det=det*A(n,n);
if abs(A(n,n))<1e-10
index=0;return;
end
for k=n:-1:1
for j=k+1:n
b(k)=b(k)-A(k,j)*x(j);
end
x(k)=b(k)/A(k,k);
end
2、输入方程组
>> A=[2/5 1/6 3/5;2/5 2/6 1/5;1/5 3/6 1/5]
A =
0.4000 0.1667 0.6000
0.4000 0.3333 0.2000
0.2000 0.5000 0.2000
>> B=[6 6 7.5]'
B =
6.0000
6.0000
7.5000
3、执行程序求解
>> [x,det,index]=Gauss(A,B)
x =
2.5000
12.0000
5.0000
det =
0.0600
index =
1
故得出结果：A、B、C三种果篮分别重2.5、12、5千克。

三、迭代法求解
1、输入Gauss-Seidel迭代法程序，保存为Gau_Seid.m。

function [x,k,index]=Gau_Seid(A,b,ep,it_max) if nargin <4 it_max=100;end
if nargin <3 ep=1e-5;end
n=length(A);k=0;
x=zeros(n,1);y=zeros(n,1);index=1;
while 1
y=x;
for i=1:n
z=b(i);
for j=1:n
if j~=i
z=z-A(i,j)*x(j);
end
end
if abs(A(i,i))<1e-5 | k==it_max
index=0;return;
end
z=z/A(i,i);x(i)=z
end
if norm(y-x,inf)<ep
break;
end
k=k+1
end
2、输入方程
（在Gauss法中已输入，故跳过）
3、执行程序求解
>> [x,k,index]=Gau_Seid(A,B,1e-5,50)
x =
2.5000
12.0000
5.0000
k =
28
index =
1
所以，A、B、C三种果篮分别为2.5、12、5千克。

四、体会
Gauss法具有中学数学解方程组的直观性，比较好理解；Gauss-Seidel迭代法是Jacobi 迭代法的改进形式，其核心思想是当新的分量求出后便立刻参与运算，可以减少计算量。

但是在使用中容易因矩阵不收敛（普半径>1）而不能使用此方法。

题目三：
一、插值法求解
1、输入分段线形插值法程序，保存为line_int_wise.m
function yi=line_int_wise(x,y,xi)
n=length(x);m=length(y);
if n~=m
error('not equal');
return;
end
for k=1:n-1
if abs(x(k)-x(k+1))<eps
error('the data is error');
return;
end
if x(k)<=xi & xi<=x(k+1)
yi=(xi-x(k+1))/(x(k)-x(k+1))*y(k)+(xi-x(k))/(x(k+1)-x(k))*y(k+1);
return;
end
end
2、输入数据
>> T=[700 720 740 760 780]
T =
700 720 740 760 780
>> v=[0.1058 0.1280 0.1462 0.1603 0.1703]
v =
0.1058 0.1280 0.1462 0.1603 0.1703
3、执行程序，求750度时的volume值
>> yi=line_int_wise(T,v,750)
yi =
0.1533
二、用最小二乘法拟合曲线
1、输入数据
(插值法中已完成，跳过)
2、绘制散点图，观察曲线形状
观察图片可知，采用一次曲线拟合即可
3、输入最小二乘法程序并构建基函数，分别保存为squar_least.m和phi_k.m function S=squar_least(x,y,n)
global i;global j;
if nargin<3 n=1; end
Phi2=zeros(n+1);
for i=0:n
for j=0:n
Phi2(i+1,j+1)=sum((phi_k(x,i)).*phi_k(x,j));
end
end
PhiF=zeros(n+1,1);
for i=0:n
PhiF(i+1)=sum((phi_k(x,i)).*y);
end
S=Phi2\PhiF;
function y=phi_k(x,k)
if k==0
y=ones(size(x));
else
y=x.^k
end
4、执行程序，获得各次项系数
>> S=squar_least(T,v,1)
S =
-0.4547
0.0008
可知：v=-0.4547+0.0008T
当T=750时，
>> v=-0.4547+0.0008.*750
v =
0.1453
三、体会
插值法容易将干扰项带入曲线，从而影响求解结果。

最小二乘法则消除了这个干扰，故在拟合时比插值法更优越。