Линейная Аппроксимация

Линейная аппроксимация

При анализе экспериментальных данных достаточно часто возникает необходимость в аппроксимации этих данных, т.е. нахождения приблизительной аналитической функции с необходимой точностью описывающей имеющиеся значения.
В ниже следующих программах на различных языках приводятся программы аппроксимации исходных данных линейной функцией по методу наименьших квадратов.

Далее приводится вариант программы на языке С (компилятор BORLAND C++ 3.1)

Show »

//***********************************************************
// Программа линейной аппроксимации табличной функции по МНК.
//***********************************************************
// Вводятся последовательно:
// - кол-во точек исходной функции
// - затем пары значений (абсцисса и ордината табличного значения)
// Внимание!!! Абсциссы не обязательно равноотстоящие.
// Контрольная точка №1: при вводе
// - 9
// - -26 66.7 -22 71 -16 76.3 -11 80.6 -5 85.7 3 92.9 10 99.4 25 113.6 42 125.1
// получаю
// 67.508 70.990 76.214 80.567 85.791 92.756 98.851 111.910 126.711
// Контрольная точка №2: при вводе
// - 3
// - 1 3 5 6 7 8
// получаю
// 2.929 6.214 7.857
// при этом а0=2.107143 а1=0.821428

#include <stdio.h>
#include <conio.h>
int n; // количество исходных точек аппроксимации
float mx[20]; // массив исходных значений абсцисс
float my[20]; // аналогично для ординат исходной табличной функции
int i; // переменная цикла
float a0,a1; // коэффициенты аппроксимирующей прямой Y=a0+a1*X
float sx=0, sy=0, sxy=0, sx2=0; // сюда буду накапливать суммы;
void main(void)
{ // 0
 clrscr(); //начальная очистка экрана
 puts("Введи количество исходных интерполируемых точек:");
 scanf("%d",&n);
 puts("Введи попарно сами исходные абсциссы (от меньших к большим) и ординаты:");
 // Ввод самих исходных значений
 for(i=0;i<=n-1; scanf("%f%f",&mx[i++],&my[i]));
 // Собственно получение искомых коээфициентов аппроксимирующей прямой
 for(i=0;i<=n-1;sx+=mx[i],sy+=my[i],sxy+=mx[i]*my[i],sx2+=mx[i]*mx[i],i++);
 a1=(sxy-sy*sx/n)/(sx2-sx*sx/n);
 a0=(sy-a1*sx)/n;
 puts("Результаты работы программы:");
 puts("номер абсцисса исх.ордината выч.ордината");
 for(i=0;i<=n-1; printf("%3d%10.3f%12.3f%14.3fn",i,mx[i++],my[i],a0+a1*mx[i]));
} // -0

Далее приводится текст программы на языке PASCAL. Проверено на компиляторах
PASCALABC 3.1 и PASCALABC.NET 2.2.

Show »

// AINTE002.CPP
//***********************************************************
// Программа линейной аппроксимации табличной функции по МНК.
//***********************************************************
// Вводятся последовательно:
// - кол-во точек исходной функции
// - затем пары значений (абсцисса и ордината табличного значения)
// Внимание!!! Абсциссы не обязательно равноотстоящие.
// Контрольная точка №1: при вводе
// - 9
// - -26 66.7 -22 71 -16 76.3 -11 80.6 -5 85.7 3 92.9 10 99.4 25 113.6 42 125.1
// получаю
// 67.508 70.990 76.214 80.567 85.791 92.756 98.851 111.910 126.711
// Контрольная точка №2: при вводе
// - 3
// - 1 3 5 6 7 8
// получаю
// 2.929 6.214 7.857
// при этом а0=2.107143 а1=0.821428

uses crt;
var 
n:integer; // количество исходных точек аппроксимации
mx: array[0..20] of real; // массив исходных значений абсцисс
my: array[0..20] of real; // аналогично для ординат исходной табличной функции
i:integer; // переменная цикла
a0,a1:real; // коэффициенты аппроксимирующей прямой Y=a0+a1*X
sx, sy, sxy, sx2:real; // сюда буду накапливать суммы;
begin // 0
 writeln('Введи количество исходных интерполируемых точек:');
 readln(n);
 writeln('Введи попарно сами исходные абсциссы (от меньших к большим) и ординаты:');
 // Ввод самих исходных значений
 // и получение искомых коээфициентов аппроксимирующей прямой
 for i:=0 to n-1 do
 begin
 read(mx[i],my[i]);
 sx:=sx+mx[i];
 sy:=sy+my[i];
 sxy:=sxy+mx[i]*my[i];
 sx2:=sx2+mx[i]*mx[i];
 end;
 a1:=(sxy-sy*sx/n)/(sx2-sx*sx/n);
 a0:=(sy-a1*sx)/n;
 writeln('Результаты работы программы:');
 writeln('номер абсцисса исх.ордината выч.ордината');
 for i:=0 to n-1 do
 writeln(i:3,mx[i]:10:3,my[i]:12:3,a0+a1*mx[i]:14:3);
end. // -0

Далее приводится текс программы на языке QBASIC в компиляторе
MS-DOS QBASIC 1.0.

Show »

'***********************************************************
' Программа линейной аппроксимации табличной функции по МНК.
'***********************************************************
' Вводятся последовательно:
' - кол-во точек исходной функции
' - затем пары значений (абсцисса и ордината табличного значения)
' Внимание!!! Абсциссы не обязательно равноотстоящие.
' Контрольная точка №1: при вводе
' - 9
' - -26 66.7 -22 71 -16 76.3 -11 80.6 -5 85.7 3 92.9 10 99.4 25 113.6 42 125.1
' получаю
' 67.508 70.990 76.214 80.567 85.791 92.756 98.851 111.910 126.711
' Контрольная точка №2: при вводе
' - 3
' - 1 3 5 6 7 8
' получаю
' 2.929 6.214 7.857
' при этом а0=2.107143 а1=0.821428

DIM n AS INTEGER ' количество исходных точек аппроксимации
DIM i AS INTEGER ' переменная цикла
DIM a0 AS SINGLE, a1 AS SINGLE' коэффициенты аппроксимирующей прямой Y=a0+a1*X
DIM sx AS SINGLE, sy AS SINGLE, sxy AS SINGLE, sx2 AS SINGLE' сюда буду накапливать суммы
INPUT "Введи количество исходных интерполируемых точек:", n
DIM mx(0 TO n - 1) AS SINGLE' массив исходных значений абсцисс
DIM my(0 TO n - 1) AS SINGLE' аналогично для ординат исходной табличной функции
PRINT "Введи попарно сами исходные абсциссы (от меньших к большим) и ординаты:"
' Ввод самих исходных значений
' и получение искомых коээфициентов аппроксимирующей прямой
FOR i = 0 TO n - 1
 INPUT mx(i), my(i)
 sx = sx + mx(i)
 sy = sy + my(i)
 sxy = sxy + mx(i) * my(i)
 sx2 = sx2 + mx(i) * mx(i)
NEXT i
a1 = (sxy - sy * sx / n) / (sx2 - sx * sx / n)
a0 = (sy - a1 * sx) / n
PRINT "Результаты работы программы:"
PRINT "номер абсцисса исх.ордината выч.ордината"
FOR i = 0 TO n - 1
 PRINT USING "### #######.### ######.### #########.###";i,mx(i),my(i),a0 + a1 * mx(i)
NEXT i
END

Далее приводится текст этой же программы, написанный на C# 4.0 в компиляторе Visual C# 2010 Express.

Show »

/***********************************************************
 Программа линейной аппроксимации табличной функции по МНК.
 Данные вводятся каждая на своей строке.
 Работа программы выглядит следующим образом:
Введи количество исходных интерполируемых точек:
3
Введи попарно сами исходные абсциссы (от меньших к большим) и ординаты каждое на
 отдельнй строке:'
1
3
5
6
7
8 
Результаты работы программы:
номер абсцисса исх.ордината выч.ордината
 001 1,000 3,000 2,929
 002 5,000 6,000 6,214
 003 7,000 8,000 7,857
***********************************************************/
using System;

namespace ConsoleApplication1
{
 class Program
 {
 static void Main(string[] args)
 {
 int n,i;
 float sx=0, sy=0, sxy=0, sx2=0, a1, a0;
 Console.WriteLine("Введи количество исходных интерполируемых точек:");
 n = int.Parse(Console.ReadLine());
 Console.WriteLine("Введи попарно сами исходные абсциссы (от меньших к большим) и ординаты:'");
 // заполнение с клавы массивов исходных данных
 float[] mx = new float[n];
 float[] my = new float[n];
 // собственно заполнение массивов в цикле необходимые вычисления
 for (i = 0; i < n; i++)
 {
 mx[i] = float.Parse(Console.ReadLine());
 my[i] = float.Parse(Console.ReadLine());
 sx+=mx[i];
 sy+=my[i];
 sxy+=mx[i] * my[i];
 sx2+=mx[i] * mx[i];
 }
 a1 = (sxy - sy * sx / n) / (sx2 - sx * sx / n);
 a0 = (sy - a1 * sx) / n;
 Console.Write("Результаты работы программы:n");
 Console.Write("номер абсцисса исх.ордината выч.ординатаn");
 for (i = 0; i < n; i++)
 Console.WriteLine(" {0:d3} {1:f3} {2:f3} {3:f3}", i+1, mx[i], my[i], a0 + a1 * mx[i]);
 Console.ReadLine(); // держит экран
 }
 }
}


Ниже приводится текст этой же программы составленный для интерпретатора Python 3.4.3.

Show »

#***********************************************************
# Программа линейной аппроксимации табличной функции по МНК.
#***********************************************************
# Вводятся последовательно:
# - кол-во точек исходной функции
# - затем пары значений (абсцисса и ордината табличного значения)
# Внимание!!! Абсциссы не обязательно равноотстоящие.
# Контрольная точка №2: при вводе
# - 3
# - 1 3 
# - 5 6
# - 7 8
# получаю
# 2.929 6.214 7.857
# при этом а0=2.107143 а1=0.821428
n= int(input('Введи кол-во исходных точек '))
print('Введи попарно сами исходные абсциссы (от меньших к большим) и ординаты:')
mmxy = [] # объявляю пустой список (в данном случае одномерный массив)
# цикл по вводу с клавы пар координат
sx=sy=sxy=sx2=0
for i in range(n): 
 # собствнно ввод с клавы через пробел
 (a,b)= list(map(float, input().split())) 
 mmxy.append(a)
 mmxy.append(b)
 sx+=a
 sy+=b
 sxy+=a*b
 sx2+=a*a
a1=(sxy - sy * sx / n) / (sx2 - sx * sx / n)
a0=(sy - a1 * sx) / n
print('Результаты работы программы:')
print(' номер абсцисса исх.ордината выч.ордината')
for i in range(n):
 print('%3d%10.3f%12.3f%14.3f'%(i+1,mmxy[2*i],mmxy[2*i+1],a0+a1*mmxy[2*i]),end='n')
input() # команда ДЕРЖИТ экран


Ниже приводится текст этой же программы составленный для компилятора Compaq Visual Fortran 2000.

Show »

! AINTE002.CPP
!***********************************************************
! Программа линейной аппроксимации табличной функции по МНК.
!***********************************************************
! Вводятся последовательно:
! - кол-во точек исходной функции
! - затем пары значений (абсцисса и ордината табличного значения)
! Внимание!!! Абсциссы не обязательно равноотстоящие.
! Контрольная точка №1: при вводе
! - 9
! - -26 66.7 -22 71 -16 76.3 -11 80.6 -5 85.7 3 92.9 10 99.4 25 113.6 42 125.1
! получаю
! 67.508 70.990 76.214 80.567 85.791 92.756 98.851 111.910 126.711
! Контрольная точка №2: при вводе
! - 3
! - 1 3 
! - 5 6
! - 7 8
! получаю
!nomer abscissa isx.ordinata wich.ordinata
! 001 1,000 3,000 2,929
! 002 5,000 6,000 6,214
! 003 7,000 8,000 7,857
! при этом а0=2.107143 а1=0.821428
real, dimension(0:20)::mmxy
integer i,n
real sx,sy,sxy,sx2! сюда буду накапливать суммы
real a0,a1 !вычисляемые коэффициенты
real a,b,kk
print*,'Wwedite kol-wo par ischodhix tochek'
read*,n
print*,'Wwedi poparno sami abscissi i ordinati:'
! Ввод самих исходных значений
! и получение искомых коээфициентов аппроксимирующей прямой
do 1 i=0,n-1
 read*,a,b
 mmxy(2*i)=a
 mmxy(2*i+1)=b
 sx = sx + a
 sy = sy + b
 sxy = sxy + a*b
 sx2 = sx2 + a*a
1 continue
a1 = (sxy - sy * sx / n) / (sx2 - sx * sx / n)
a0 = (sy - a1 * sx) / n
print*,'nomer abscissa isx.ordinata wich.ordinata'
do 2 i = 0,n - 1
 kk=a0 + a1 * mmxy(2*i)
 write(*,7)i,mmxy(2*i),mmxy(2*i+1),kk
7 format(' ',I3,F10.3,F12.3,F14.3)
2 continue 
end