MATLAB-绘图，用高斯消元法求解线性方程组

作业4-第三周周六

第1，2题见作业3

第3题

绘制如下方程给出的椭圆：
(x-2)²+2*(y-x-1)² =3，
把y看做x的函数解方程。这样做将得到两个方程，一个是椭圆的上半部分，另一个是椭圆的下半部分。绘制这两个关于x的方程，得到椭圆。

clear
syms x y;
S=solve((x-2)^2+2*(y-x-1)^2==3,y)

S =
 
 x - (2^(1/2)*(- x^2 + 4*x - 1)^(1/2))/2 + 1
 x + (2^(1/2)*(- x^2 + 4*x - 1)^(1/2))/2 + 1

clear
x=0:0.01:4;
f1=x-(2^(1/2)*(- x.^2 + 4.*x - 1).^(1/2))/2 + 1;
f2=x+(2^(1/2)*(- x.^2 + 4.*x - 1).^(1/2))/2 + 1;
plot(x,[f1;f2]);

---

3.1.解方程

代码：

clear
syms x y;
S=solve((x-2)^2+2*(y-x-1)^2==3,y)

结果：

S =
 
 x - (2^(1/2)*(- x^2 + 4*x - 1)^(1/2))/2 + 1
 x + (2^(1/2)*(- x^2 + 4*x - 1)^(1/2))/2 + 1

3.2.求x的定义域

代码：

clear
syms x
S=solve((- x^2 + 4*x - 1)^(1/2)==0,x)

结果：

S =
2 - 3^(1/2)
3^(1/2) + 2

3.3. 绘图

代码：

clear
x=2 - 3^(1/2):0.001: 3^(1/2) + 2
f1=x-(2^(1/2)*(- x.^2 + 4.*x - 1).^(1/2))/2 + 1;
f2=x+(2^(1/2)*(- x.^2 + 4.*x - 1).^(1/2))/2 + 1;
plot(x,[f1;f2]);

结果：

第4题

将两个设计参数x和y定为0<x<5和0<y<5。产品的成本为：

用mesh命令绘图，近似地找出使成本最低的最优参数，并写出最低成本值。

代码：

x=linspace(0,5,500);
y=linspace(0,5,500);
[X,Y]=meshgrid(x,y);
f=X.^2-8.*X+Y.^2-6.*Y-0.1.*X.*Y+50;
mesh(X,Y,f);
min_f=min(f,[],'all')
[m,n]=find(f==min_f)
x0=5/500*m
y0=5/500*n

结果：

min_f =
   23.7343
m =
   321
n =
   416
x0 =
    3.2100
y0 =
    4.1600

所以，使成本最低的最优参数：x0 =3.2100，y0 =4.1600；最低成本值：min_f =23.7343。

5.用contour命令重做上题。

x=linspace(0,5,500);
y=linspace(0,5,500);
[X,Y]=meshgrid(x,y);
f=X.^2-8.*X+Y.^2-6.*Y-0.1.*X.*Y+50;
contour3(X,Y,f,500)
min_f=min(f,[],'all')
[m,n]=find(f==min_f)
x0=5/500*m
y0=5/500*n

结果：

min_f =
   23.7343
m =
   321
n =
   416
x0 =
    3.2100
y0 =
    4.1600

所以，使成本最低的最优参数：x0 =3.2100，y0 =4.1600；最低成本值：min_f =23.7343。

---

6.编写消去法解线性方程组的程序，并求解

0.002x1+2x2+2x3 =0.4
x1+0.78125x2 =1.3816
3.996x1+5.5625x2+4x3=7.4178

6-1.编写的高斯消元函数：

代码：

function x = Gauss(Matrix,n,b)
% 高斯消元法
% 输入系数矩阵，系数矩阵的维数和常向量

% Step1 顺序消元

% Step1-1 排除无解或有无数个解的情况
if det(Matrix)==0
    error('奇异矩阵')
    return
end

% Step1-2 消元过程
for i=1:n-1
% Step1-2-1 当遇到矩阵的消元因子等于0的情况时，进行“行交换”
    if Matrix(i,i)==0 % 如果矩阵的消元因子等于0
        for m=i+1:n %查找下方第一个非零元素
            if Matrix(m,i)~=0 
                for p=i:n  %从第i列开始，交换第m行和第i行的元素
                    tmp=Matrix(i,p);
                    Matrix(i,p)=Matrix(m,p);
                    Matrix(m,p)=tmp;
                end
                t=b(i);b(i)=b(m);b(m)=t; % 交换向量b第m行和第i行的元素
                break % 跳出查找非零元素的循环
            end
        end
end

% Step1-2-2 对第i行下面的每一行进行初等行变换，最终使得（i，i）下方的元素都为0
    for j=i+1:n 
        if Matrix(j,i)==0 % 如果M(j,i)已经为0，则该行不需要进行行变换，继续下一行（这个步骤没有也可以）
            continue
        end
        fac1=Matrix(j,i)/Matrix(i,i); % 求第i行元素需要乘的系数
        for k=i:n % 对第j行的每个元素进行初等行变换
            Matrix(j,k)=Matrix(j,k)-Matrix(i,k)*fac1;
        end
        b(j)=b(j)-b(i)*fac1; % 对向量b的第j行进行同样的处理
    end
end

% Step2 回代求解
x(n)=b(n)/Matrix(n,n); % 求出x的最后一行的元素
for i=n:-1:2
    for j=i-1:-1:1  % 消去Matrix(i,i)正上方的元素
        fac2=Matrix(j,i)/Matrix(i,i);
        Matrix(j,i)=0;
        b(j)=b(j)-b(i)*fac2; % 对b(j)进行相应的处理
    end
    x(i-1)=b(i-1)/Matrix(i-1,i-1); % 求出x(i-1)
end

end

6-2. 求解方程

代码：

A=[-0.002 2 2;1 0.78125 0;3.996 5.5625 4];
b=[0.4 1.3816 7.4178]';
x=Gauss(A,3,b)
format long

结果：

x =

   1.927300000000076  -0.698495999999967   0.900423299999967

第十题

A=zeros(20);
for i=1:20
    A(i,i)=5
end
for i=1:19
    A(i,i+1)=-2;
    A(i+1,i)=-2;
end
for i=1:18
    A(i,i+2)=1;
    A(i+2,i)=1
end

b=[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
x=Gauss(A,20,b)
format long

结果：

x =

  列 1 至 5

   0.242121373548085   0.094391354624679  -0.021824158491067  -0.031362343009361  -0.006942557862112

  列 6 至 10

   0.004886927715767   0.003586117214438   0.000214823135181  -0.000784526508185  -0.000357861979163

  列 11 至 15

   0.000050437638520   0.000106309021306   0.000029232399870  -0.000014343050585  -0.000012667696430

  列 16 至 20

  -0.000001464361499   0.000002491258113   0.000001311981447  -0.000000093354238  -0.000000299737984