MATLAB 与C语言实现追赶法 （托马斯算法）

算法的原理就是类似于LU分解，然后求解线性方程组。具体可以参考https://www.cnblogs.com/xpvincent/archive/2013/01/25/2877411.html

方法一

a = [0,1,1,1]';
b = [4,4,4,4]';
c = [1,1,1]';
x = [1,1,1,0]';
m = zeros(3,1);
n = 4;
c(1) = c(1)/b(1);
x(1) = x(1)/b(1);
for i = 2:n-1
     c(i) = c(i)/(b(i)-c(i-1)*a(i));
end
 for i = 2:n
    x(i) = (x(i)-x(i-1)*a(i))/(b(i)-c(i-1)*a(i));   
end

for i = n-1:-1:1
    x(i) = x(i)-c(i)*x(i+1);
end

方法二 

clear all;clc;
prompt = 'What is the original value? ';
n = input(prompt);
n
a=zeros(1,n);b=zeros(1,n);c=zeros(1,n-1);
A=zeros(n,n);
for i=1:n
    for j=1:n
        if j==i
            A(i,j)=4;
        elseif (j-i==-1)||(j-i==1)
            A(i,j)=1;
        else
            A(i,j)=0;
        end
    end
end
fprintf('系数矩阵：');A
fprintf('主对角线元素：');b=diag(A)'
fprintf('第1条对角线元素：');c=diag(A,1)'
fprintf('第-1条对角线元素：');a(2:n)=c
d=[1,1,1,0];%zeros(1,n);d(1)=1;d(n)=((-1)^(n+1));
fprintf('给定的系数矩阵(转置后)：');d
y=zeros(n,1);x=y;
u=zeros(1,n);l=u;%u矩阵U的主对角线元素，l矩阵L的第-1条对角线元素
%---------追过程--------
u(1)=b(1);y(1)=d(1);
for i=2:n
    l(i)=a(i)/u(i-1);
    u(i)=b(i)-l(i)*c(i-1);
    y(i)=d(i)-l(i)*y(i-1);
end
fprintf('中间解向量：');y
fprintf('单位下三角阵L的第-1条对角线元素：');l 
fprintf('上三角阵U的主对角线元素：')
U=zeros(n);
L=eye(n);
for i=1:n-1
    L(i+1,i)=l(i+1);
end
fprintf('单位下三角阵：');%单位下三角阵
 
for i=1:n-1
    U(i,i)=u(i);
    U(i,i+1)=c(i);
end
U(n,n)=u(n);
fprintf('上三角阵：'); 
%---------赶过程--------
x(n)=y(n)/u(n);
for i=n-1:-1:1
    x(i)=(y(i)-c(i)*x(i+1))/u(i);
end
fprintf('解向量：');x

C代码

#include<stdio.h>

void tdma(float x[], const size_t N, const float a[], const float b[], float c[]);
int main()
{
    float x[4] = { 1,1,1,0 };
    int i;
    const size_t N = 4;
    const float a[4] = {0,1,1,1};
    const float b[4] = { 4,4,4,4 };;
    float c[3] = { 1,1,1 };

    tdma(x ,N,a ,b ,c );
    for ( i = 0; i < N; i++)
    {
        printf("%f ", x[i]);
    }  
    return 0;
} 

void tdma(float x[], const size_t N, const float a[], const float b[], float c[])
{
    size_t n;

    c[0] = c[0] / b[0];
    x[0] = x[0] / b[0];

    for (n = 1; n < N-1; n++) {
         
        c[n] = c[n] / (b[n] - a[n] * c[n - 1]);
    }
    for (n = 1; n < N ; n++) {
         
        x[n] = (x[n] - a[n] * x[n - 1]) /(b[n] - a[n] * c[n - 1]);
    }

    for (n = N -2; n-- > 0; )
        x[n] = x[n] - c[n] * x[n + 1];
}