Линейная аппроксимация
При анализе экспериментальных данных достаточно часто возникает необходимость в аппроксимации этих данных, т.е. нахождения приблизительной аналитической функции с необходимой точностью описывающей имеющиеся значения.
В ниже следующих программах на различных языках приводятся программы аппроксимации исходных данных линейной функцией по методу наименьших квадратов.
Далее приводится вариант программы на языке С (компилятор BORLAND C++ 3.1)
//*********************************************************** // Программа линейной аппроксимации табличной функции по МНК. //*********************************************************** // Вводятся последовательно: // - кол-во точек исходной функции // - затем пары значений (абсцисса и ордината табличного значения) // Внимание!!! Абсциссы не обязательно равноотстоящие. // Контрольная точка №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.
// 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.
'*********************************************************** ' Программа линейной аппроксимации табличной функции по МНК. '*********************************************************** ' Вводятся последовательно: ' - кол-во точек исходной функции ' - затем пары значений (абсцисса и ордината табличного значения) ' Внимание!!! Абсциссы не обязательно равноотстоящие. ' Контрольная точка №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.
/*********************************************************** Программа линейной аппроксимации табличной функции по МНК. Данные вводятся каждая на своей строке. Работа программы выглядит следующим образом: Введи количество исходных интерполируемых точек: 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.
#*********************************************************** # Программа линейной аппроксимации табличной функции по МНК. #*********************************************************** # Вводятся последовательно: # - кол-во точек исходной функции # - затем пары значений (абсцисса и ордината табличного значения) # Внимание!!! Абсциссы не обязательно равноотстоящие. # Контрольная точка №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.
! 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