北邮Matlab选修 Week 7 线性方程组求解

从课程所讲授的线性方程组求解方法中任选一种方法实现对下面的线性方程组的求解

问题描述

解方程

程序设计

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12


%solveByCramer.m
function x = solveByCramer(A, b)
    n = length(b);
    x = zeros(n, 1);
    detA = det(A);

    for i = 1:n
        Ai = A;
        Ai(:,i) = b;
        x(i) = det(Ai) / detA;
    end
end

计算结果

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26


>>  A=[1,2,3,4;4,3,2,1;1,3,2,4;4,1,3,2]

A =

     1     2     3     4
     4     3     2     1
     1     3     2     4
     4     1     3     2

>>  b=[5;4;3;2]

b =

     5
     4
     3
     2

>>  result=solveByCramer(A,b)

result =

   -1.8000
    1.8667
    3.8667
   -2.1333

电路问题

问题描述

一种大型输电网络可以简化为下图所示电路，其中$R_i$表示负载电阻，$r_i$表示线路内阻，设电源电压为$V$。若$R_i=6$，$r_i=1$，$V=18$，求出各个负载上的电流$I_1,I_2,I_3, \cdots , I_n$及总电流$I_0$。

数学模型

物理关系

对于支路电源和$R_1$，$r_1$，有： $$ V=U_{R_1}+U_{r_1} $$

即： $$ V=I_1 \times R_1 + I_0 \times r_1 $$

可化为： $$ I_0+6 I_1 = 18 $$ (共 1 个方程)
对于支路$R_{k+1}$，$r_{k+1}$，$R_k$，有： $$ U_{R_k}=U_{R_{k+1}}+U_{r_{k+1}} (1 \leq k < n) $$

即： $$ I_k \times R_k = I_{k+1} \times R_{k+1} + \sum_{j=k+1}^{n} I_j \times r_k $$

可化为： $$ -6I_k+7I_{k+1}+\sum_{j=k+2}^{n} I_j = 0 $$ (共 n-1 个方程)
根据电流关系，有： $$ I_0=\sum_{j=1}^{n} I_j $$ (共 1 个方程)

以上共计n+1个方程，对应$I_0$到$I_n$n+1个未知数。

转化成数学问题

转化成矩阵运算：
解以上方程组可得解向量$I_0$~$I_n$

程序设计

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31


%Circuit.m
function result = Circuit(n)
    % -6*x_{k} + 7*x_{k+1} + (x_{k+2} + x_{k+3} + ... + x_{n}) = 0
    A = zeros(n+1, n+1);
    b = zeros(n+1, 1);
    b(n+1) = 18;
    % Fill in the matrix A according to the recursive relation
    for k = 2:n
        A(k, k) = -6;
        A(k, k+1) = 7;
        if k+2 <= n+1
            A(k, k+2:n+1) = 1;
        end
    end
    A(1,1)=1;
    A(1,2:n+1)=-1;
    A(n+1,1)=1;
    A(n+1,2)=6;
    result = solveByCramer(A,b);
end
function x = solveByCramer(A, b)
    n = length(b);
    x = zeros(n, 1);
    detA = det(A);

    for i = 1:n
        Ai = A;
        Ai(:,i) = b;
        x(i) = det(Ai) / detA;
    end
end

计算结果

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15


>> Circuit(10)

ans =

    5.9970
    2.0005
    1.3344
    0.8907
    0.5955
    0.3995
    0.2702
    0.1858
    0.1324
    0.1011
    0.0867

最后修改于 2024-08-05