4392

Численное решение уравнений в языке С++

Контрольная

Информатика, кибернетика и программирование

Численное решение уравнений в языке С++ Теоретические основы Предположим, нам нужно решить кубическое уравнение Это означает, что нужно найти корни уравнения – такие числа, которые обращают уравнение в ноль...

Русский

2012-11-18

168.5 KB

70 чел.

Численное решение уравнений в языке С++

  1.  Теоретические основы

Предположим, нам нужно решить кубическое уравнение

.     (8.1)

Это означает, что нужно найти корни уравнения – такие числа, которые обращают уравнение в ноль, если приравнять неизвестную переменную x какому-либо из этих чисел. Основная теорема алгебры утверждает, что корней может быть несколько, причем, некоторые из них могут являться комплексными числами. Обычно практический интерес представляют только действительные корни. Убедиться в наличии действительных корней в уравнении (8.1) можно путем построения графика функции

.    (8.2)

Такой график приведен на рис. 8.1.

Рис. 8.1. График функции

Как можно видеть, действительный корень уравнения (8.1) располагается где-то в диапазоне между 0,6 и 0,8. Если известно, что в заданном диапазоне находится только один корень уравнения, то функции от крайних точек диапазона, как правило, имеют разные знаки. Например, в случае функции, изображенной на графике рис. 8.1, имеет положительное значение, а  – отрицательное значение. Можно утверждать и обратное: если значения функции на краях интервала имеют разные знаки, то внутри интервала содержится корень уравнения. Однако если функция имеет разрывы в некоторых точках, это правило не справедливо. Например, функция  изменяет знак в точках , , как показано на рис. 8.2. Однако в этих точках нет корней уравнения, поскольку функция  не пересекает ось . Такие точки называются сингулярными.

Рис. 8.2. График функции

  1.  Метод простого перебора

Самым простым методом решения алгебраических уравнений является метод перебора в заданном интервале. Программа, которая реализует данный метод, представлена в листинге 8.1. Помимо основной функции main() она содержит также функции RootSearch()и Func(). В функции main() пользователь вводит исходные данные, после чего вызывается функция RootSearch(), которая непосредственно и вычисляет корень уравнения.

Листинг 8.1. Решение уравнения методом перебора

#include <iostream.h>

/* Ниже следуют прототипы функций */

void RootSearch ();

double Func (double);

int n;

double p[10];

double a, b, dx;  // Глобальные переменные

void main ()

{

cout<<"Input degree of polynomial: ";

cin>>n;

cout<<endl;

cout<<"Input coefficients of function:\n";

cout<<endl;

for (int i=n; i>=0; i--) cin>>p[i];

cout<<endl;

cout<<"Input left limit: ";

cin>>a;

cout<<endl;

cout<<"Input right limit: ";

cin>>b;

cout<<endl;

cout<<"Input step: ";

cin>>dx;

cout<<endl;

RootSearch();

char Res;

cin>>Res;

}

void RootSearch ()

{

double x1, x2, f1, f2; //Локальные переменные

 x1=a;

x2=a+dx;

f1=Func(x1);

f2=Func(x2);

while (f1*f2>0.0)

{

 if (x1>=b)

 {

  cout<<"Root is not bracketed in

("<<a<<","<<b<<")"<<"\n";

  return;

 }

 x1=x2;

 f1=f2;

 x2=x1+dx;

 f2=Func(x2);

}

cout<<"x1: "<<x1<<"\n";

cout<<"x2: "<<x2<<"\n";

}

double Func (double x)

{

double Sum=0;

double y=1;

for (int i=0; i<(n+1); i++)

{

 Sum=Sum+p[i]*y;

 y=y*x;

 }

return Sum;

}

РЕЗУЛЬТАТ

Функция Func() используется функцией RootSearch() как подпрограмма. Во второй и третьей строчках программы заявлены прототипы функций. Далее объявляются глобальные переменные, которые доступны для всех функций программы. В языке С++ глобальные переменные почти никогда не используются, поскольку существует возможность, что какая-либо из функций может изменить их значение. Альтернативой глобальным переменным являются статические переменные-члены, которые доступны всем экземплярам своего класса (термин из области объектно-ориентированного программирования). Они представляют собой компромисс между глобальными данными, доступными всем элементам программы, и данными-членами, доступными только объектам этого класса. В данной программе глобальные переменные используются в качестве примера. Угроза их несанкционированного изменения отсутствует. В отличие от глобальных переменных локальные переменные объявляются внутри функций. После выполнения функцией своей задачи ее локальные переменные стираются из памяти.

  1.  Метод половинного деления

Метод половинного деления называется также методом дихотомии (от греческого – разделяю на две части). Англоязычные народы называют его method of bisection. Этот метод не самый скорый, однако наиболее надежный. Он использует следующий принцип: если корень находится в интервале , то произведение .

При использовании метода дихотомии вначале определяют точку в середине отрезка: . Затем находят значение функции в этой точке . Если , то корень заключен в интервале , а если , то – в интервале . Далее операцию деления повторяют с тем отрезком, в котором может находиться корень. Итерации продолжают до тех пор, пока интервал поиска не станет меньше заданной малой величины :

.

Исходный интервал  уменьшается до  после одного деления. А после  делений уменьшается до величины . Приравнивая , и решая это уравнение относительно , получим

.     (8.3)

В листинге 8.2 представлена программа, решающая заданное уравнение методом половинного деления. Для определения количества делений используется формула (8.3). В данной программе помимо основной функции main() содержатся также функции Bisect()и Func(). В отличие от предыдущего примера (листинг 8.1) данные для вычислений вводятся в самой функции Bisect(), а функция main() используется только для вызова функции Bisect(). Такой подход считается более приемлемым.

Входной аргумент filter предназначен для |контролируфильтрации возможной сингулярности. Устанавливая filter=1, мы заставляем программу проверять, уменьшается ли величина функции  с каждым делением интервала пополам. Если нет, искомый «корень», является не корнем, а сингулярностью. Если есть уверенность, что сингулярность в данной задаче отсутствует, следует установить: filter=0.

Листинг 8.2. Решение уравнения методом дихотомии

#include <iostream.h>

#include <math.h>

int Bisect();  // Прототипы функций

double Func(double, int, double);

 

void main()

{

Bisect();

}

int Bisect()

{

int i, n, k, filter;

double p[10];

double a, b, x1, x2, x3, f1, f2, f3, tol, root;

cout<<"Input degree of polynomial: ";

cin>>n;

cout<<endl;

cout<<"Input coefficients of function:\n";

cout<<endl;

for (int i=n; i>=0; i--) cin>>p[i];

cout<<endl;

cout<<"Input left limit: ";

cin>>a;

cout<<endl;

cout<<"Input right limit: ";

cin>>b;

cout<<endl;

cout<<"Input tol: ";

cin>>tol;

cout<<endl;

//filter=singularity filter: 0=off (default), 1=on.

cout<<"Input filter (0 or 1): ";

cin>>filter;

cout<<endl;

f1=Func(p, n, a);

if (f1==0.0)

{

 root=a;

 cout<<"Root: "<<root<<"\n";

 char Res;

 cin>>Res;

 return 0;

}

f2=Func(p, n, b);

if (f2==0.0)

{

 root=b;

 cout<<"Root: "<<root<<"\n";

 char Res;

 cin>>Res;

 return 0;

}

if (f1*f2>0)

{

 cout<<"Root is not bracketed in ("<<a<<", "<<b<<")"<<"\n";

 char Res;

 cin>>Res;

 return 0;

}

k=ceil(log(fabs(b-a)/tol)/log(2.0));

x1=a;

x2=b;

for (

 i=1; i<(k+1); i++)

{

 x3=0.5*(x1+x2);

 f3=Func(p, n, x3);

 if ((filter==1)&&(fabs(f3)>fabs(f1))&&(fabs(f3)>fabs(f2)))

 {

  cout<<"Root is not bracketed in ("<<a<<",

"<<b<<")"<<"\n";

  char Res;

  cin>>Res;

  return 0;

 }

 if (f3==0.0)

 {

  root=x3;

  cout<<"Root: "<<root<<"\n";

  char Res;

  cin>>Res;

  return 0;

 }

 if (f2*f3<0.0)

 {

  x1=x3;

  f1=f3;

 }

 else

 {

  x2=x3;

  f2=f3;

 }

}

root=0.5*(x1+x2);

cout<<"Root: "<<root<<"\n";

char Res;

cin>>Res;

return 0;

}

double Func (double p[], int n, double x)

{

double Sum=0;

double y=1;

for (int i=0; i<(n+1); i++)

{

 Sum=Sum+p[i]*y;

 y=y*x;

}

return Sum;

}

  1.  Метод Ньютона-Рафсона

Алгоритм Ньютона-Рафсона является наиболее известным методом отыскания корней уравнения по причине своей простоты и высокого быстродействия. Гладкую и непрерывную функцию можно разложить в ряд Тейлора

.

Если является корнем уравнения , то можем записать

.

Приравнивая остаток  нулю, получим формулу Ньютона-Рафсона:

.     (8.4)

Эта формула позволяет находить корень уравнения (если он существует на заданном отрезке) с заданной ошибкой за конечное число шагов. Программа, решающая заданное уравнение методом Ньютона-Рафсона, представлена в листинге 8.3.

Листинг 8.3. Решение уравнения методом Ньютона-Рафсона

#include <iostream.h>

#include <math.h>

void NewtonRaphson(int, double);

double Func(int, double);

double dFunc(int, double);

double p[10];

void main (){

int n;

double x;

cout<<"Input degree of polynomial: ";

cin>>n;

cout<<endl;

cout<<"Input coefficients of function:\n";

cout<<endl;

for (int i=n; i>=0; i--) cin>>p[i];

cout<<endl;

cout<<"Input x: ";

cin>>x;

cout<<endl;

NewtonRaphson(n, x);

char Res;

cin>>Res;

}

 

void NewtonRaphson(int n, double x)

{

const double tol=0.000001;

double dx, root;

for (int i=1; i<31; i++)

{

 dx=-Func(n, x)/dFunc(n, x);

 x=x+dx;

 if (fabs(dx)<tol)

 {

  root=x;

  cout<<"Root: "<<root<<endl;

  return;

 }

}

cout<<"The root isn't found";

}

double Func (int n, double x)

{

double Sum=0;

double y=1;

for (int i=0; i<(n+1); i++)

{

 Sum=Sum+p[i]*y;

 y=y*x;

}

return Sum;

}

double dFunc (int n, double x)

{

double Sum=0;

double y=1;

for (int i=1; i<(n+1); i++)

{

 Sum=Sum+i*p[i]*y;

 y=y*x;

}

return Sum;

}

РЕЗУЛЬТАТ

 

На рис. 8.3 показан график функции

.  (8.5)

Как можно видеть, эта функция на интервале (0, 5) имеет два корня. Причем, первый корень совпадает с локальным экстремумом функции, при этом значения функции слева и справа от этого корня имеют один и тот же знак. Поэтому метод дихотомии не позволит определить этот корень. Однако метод Ньютона-Рафсона вполне способен справиться с данной проблемой.

                    

Рис. 8.3. График функции

Численные методы позволяют найти лишь один корень – даже в том случае, если уравнение имеет несколько действительных корней. Например, если в качестве исходных данных мы возьмем полином (8.5), и в качестве начального значения примем х=3, то получим корень х=2,1. Если же мы возьмем начальное значение х=3,9; то получим второй корень х=4.

  1.  Решение систем линейных уравнений методом Гаусса

Метод Гаусса относится к наиболее известным методам решения систем линейных уравнений. Он делится на две части: фазу исключения и фазу обратной подстановки. Рассмотрим этот метод на примере, решив следующую систему из трех уравнений.

,    (a)

,    (b)

.    (c)

 Фаза исключения. Эта фаза заключается в умножении одного уравнения на одно и то же число , и затем вычитания полученного произведения из другого уравнения.  В результате количество членов преобразованного уравнения уменьшается. Символически это можно записать так

.

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

,

.

После этой трансформации получим следующую систему уравнений

,    (a)

,    (b)

.    (c)

Далее в качестве стержневого уравнения выбираем уравнение (b) и преобразуем уравнение (c) по следующей схеме:

.

В результате получим

,    (a)

,    (b)

.      (c)

Фаза обратной подстановки. Поскольку мы получили систему с верхней треугольной матрицей, то вначале находим неизвестную с наибольшим индексом, а затем – остальные неизвестные:

,

  ,

.

Алгоритм метода Гаусса. Рассмотрим систему в следующем виде.

Пусть i-тая строка является текущей строкой, лежащей ниже стержневой, ее мы будем трансформировать. Для исключения  из этой строки  элемента, нужно, чтобы сомножитель  был равен отношению . Получаем следующий алгоритм:

,

    .

После приведения матрицы к треугольному виду k-я строка имеет вид:

.

Поэтому неизвестные переменные можно найти по формуле:

.

 

Листинг 8.4. Решение систем линейных уравнений методом Гаусса

#include <iostream.h>

void input();

void output();

void gauss();

double A[100][100], b[100];

int n;

void main() {

input();

gauss();

cout<<endl;

cout<<"Result:"<<endl;

output();

char res;

cin>>res;

}

void input() {

cout<<"Input matrix size: ";

cin>>n;

cout<<"Input matrix:"<<endl;

cout<<endl;

cout<<"A = ";

for (int i=1; i<=n; i++) {

 for (int j=1; j<=n; j++) {

  cin>>A[i][j];

 }

}

cout<<endl;

cout<<"Input vector:"<<endl;

cout<<endl;

cout<<"b = ";

for (i=1; i<=n; i++)

 cin>>b[i];

}

void output() {

cout<<endl;

for (int i=1; i<=n; i++) {

 for (int j=1; j<=n; j++) {

  cout<<A[i][j]<<"\t";

 }

 cout<<endl;

}

cout<<endl;

for (i=1; i<=n; i++)

 cout<<b[i]<<"\n";

}

void gauss() {

double lambda;

int i, j, k;

// Elimination phase

for (k=1; k<n; k++) {

 for (i=k+1; i<=n; i++) {

  if (A[i][k]!=0) {

   lambda=A[i][k]/A[k][k];

   b[i]=b[i]-lambda*b[k];

   for (j=k; j<=n; j++) {

    A[i][j]=A[i][j]-lambda*A[k][j];

   }

  }

 }

}

// Back substitution phase (x=b, Sum=lambda)

for (k=n; k>0; k--) {

 lambda=0;

 for (j=k+1; j<=n; j++) {

  lambda=lambda+A[k][j]*b[j];

 }

 b[k]=(b[k]-lambda)/A[k][k];

 }

}

РЕЗУЛЬТАТ


 

А также другие работы, которые могут Вас заинтересовать

23101. Хвильові властивості частинок. Хвилі де Бройля 5.02 MB
  Хвилі де Бройля. Тобто інколи відбувається прояв як хвилі інколи як частинки. Тоді можна отримати вираз для хвилі де Бройля . Оберемо напрям вздовж за напрям розповсюдження хвилі де фаза хвилі що пересувається у просторі з фазовою швидкістю що шукається з умови що переміщується так щоб фаза залишалась постійною .
23102. Принципова схема лазера. Властивості лазерного випромінювання. Типи лазерів та їх застосування 51.5 KB
  При падінні хвилі з власною частотою переходу системи: змінюються заселеності рівнів N1 i N2 кількість атомів в одиниці об’єму що знаходяться на 1 та на 2 енергетичних рівнях відповідно. dN12=BN1dt ; кількість частинок що перейшли з 1 рівня на 2 dN21= AN2dt BN2dt кількість частинок що перейшли з 2 рівня на 1 де Акоеф. Крім того в стаціонарному режимі при умові термодинамічної рівноваги виконуються рівняння: N1N2=N=const кількість частинок в системі є сталою. В дворівневій системі не можна забезпечити умову N2 N1 бо навіть в...
23103. Рівняння Шредингера. Інтерпретація хвильової функції 49 KB
  Рівняння Шредингера. Для цього необхідне рівняння: 1. Рівняння повинно бути лінійним і однорідним хвиля задовольняє принц. Це рівняння Шредингера.
23104. Співвідношення невизначеності Гейзенберга, приклади його проявів 74.5 KB
  Нехай стан частинки опивується хв. Остаточно Співвідношення невизначеностей проявляється при будьякій спробі вимірювання точного положення або точного імпульса частинки. Виявляється що уточнення положення частинки впливає на те що збільшується неточність в значенні імпульса і навпаки. Часто втрачає зміст ділення повної енегрії частинкияк квантового об’єкту на потенціальну і кінетичну .
23105. Сестринский процесс при холециститах 25.25 MB
  Воспаление желчного пузыря регистрируется почти у 10% населения планеты, причем в 3-4 раза чаще холециститом страдают женщины. Большинство людей не следят за своим рационом, ведут сидячий образ жизни.
23106. Теорія молекули водню. Обмінна взаємодія 371 KB
  Оскільки гамільтоніан не залежить від спінових змінних то хвильова функція зображається добутком спінової функції на просторову . За допомогою хвильової функції знаходимо середнє значення повного гамільтоніана системи: де кулонівський інтеграл К характаризує ел. наближені хвильові функції Кулонівський інтеґрал К є малим числом і головну роль відіграє обмінний інтеґрал який у ділянці малих є додатною величиною а при змінює знак. Таким чином для симетричної просторової функції є можливим зв'язаний стан системи і теорія...
23107. Прискорювачі заряджених частинок та принципи їх роботи 62.5 KB
  При непрямих методах прискорення електричне поле індукується змінним магнітним полем або використовується змінне електричне поле у вигляді біжучих або стоячих хвиль. Ідея прискорення заряджених частинок електричним полем яке породжується змінним магнітним полем. Основна складова – потужний електромагніт обмотка якого живиться змінним струмом з частотою сотні МГц. При зміні маг потока з’являється вихрове ел поле і на кожний електрон в камері діє сила eE.
23108. Общая характеристика экономики государственного сектора 262 KB
  Под государственным сектором экономики страны понимают сектор, представляющий и обслуживающий интересы всего населения. Государство является основным институтом, организующим и координирующим взаимоотношения граждан и социальных групп в стране и обеспечивающим условия для их совместной деятельности
23109. Сучасні уявлення про ядерні сили. Моделі атомного ядра 136.5 KB
  За сучасними поглядами сили між нуклонами є виявом сильної кваркглюонної взаємодії. Така частинканосій сильної міжкваркової взаємодії називається глюоном. При взаємодії глюонів з кварками колір кварків змінюється. Аромат кварків їхній електричний та баріонний заряди не змінюються тобто колір є найбільш важливою властивістю кварків при взаємодії.