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];

 }

}

РЕЗУЛЬТАТ


 

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

60874. Язык разметки гипертекста HTML 50 KB
  Ожидаемые результаты: В конце урока ученики смогут: Дать определение таким понятиям: Гипертекст; WEBстраница WEBсайт; WEBдизайн; Язык разметки гипертекста HTML; Использовать команды языка HTML для форматирования текста; Ориентировочный план урока Актуализация опорных знаний Как вы можете определить понятие Интернет Что такое служба WWW Как называются документы которые мы можем просматривать в программебраузере Изложение нового материала. Систематизация понятий Гипертекст; WEBстраница WEBсайт;...
60875. Нахождение площади фигур при решении практических задач 542 KB
  Четырехугольник у которого противолежащие стороны параллельны; По вертикали: Параллелограмм у которого все углы прямые; Прямоугольник у которого все стороны равны...
60878. Лексика. Росіянізм, калька, покруч 95.5 KB
  Русизми це окремі російські слова найчастіше з українізованою вимовою наприклад: тормозити замість гальмувати порівняйте рос. Затребуваний і кальковані від нього слова: запитаний запотребуваний затребований восребуваний витребуваний.
60880. Опера М.И.Глинки Иван Сусанин 25.58 KB
  Семья Медичи славилась своим пристрастием к искусству, особенно любовью к театру и музыке, и все предвкушают, что на свадьбе смогут увидеть и услышать нечто совершенно необычное. Дворец уже полон гостей. Здесь собрался цвет итальянского искусства
60881. Общие методы решения уравнений 64.5 KB
  Цель урока: формировать способность к обобщению знаний учащихся о методах решения уравнений; способствовать дальнейшей отработке знаний и умений применять общие методы при решении уравнений.