4392

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

Контрольная

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

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

Русский

2012-11-18

168.5 KB

96 чел.

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

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

 }

}

РЕЗУЛЬТАТ


 

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

74777. Сложение взаимно-перпендикулярных колебаний. Фигуры Лиссажу 70 KB
  Если частоты складываемых взаимно перпендикулярных колебаний различны, то замкнутая траектория результирующего колебания довольно сложна. Замкнутые траектории, прочерчиваемые точкой, совершающей одновременно два взаимно перпендикулярных колебания, называются фигурами Лиссажу.
74778. Скорость света. Опыт Майкельсона. Принцип относительности Эйнштейна. Скорость света 18.63 KB
  Скорость света относится к фундаментальным физическим постоянным которые характеризуют не просто отдельные тела а свойства мира в целом. По современным представлениям скорость света в вакууме предельная скорость движения частиц и распространения взаимодействий.
74779. Преобразования Лоренца. Релятивистское изменение длин и промежутков времени. “Парадокс близнецов” 44.5 KB
  Преобразованиями Лоренца в физике в частности в специальной теории относительности СТО называются преобразования которым подвергаются пространственно-временные координаты xyzt каждого события при переходе от одной инерциальной системы отсчета ИСО к другой.
74780. Релятивистский закон сложения скоростей. Изменение массы со скоростью. Связь массы и энергии 56 KB
  Произведя соответствующие преобразования получаем релятивистский закон сложения скоростей специальной теории относительности: Закон взаимосвязи массы и энергии Найдем кинетическую энергию релятивистской частицы.
74781. Термодинамический и статистический подход к изучению поведения систем. Термодинамические параметры. Статистическое и термодинамическое определение абсолютной температуры 30.5 KB
  Законы поведения огромного числа молекул, являясь статистическими закономерностями, изучаются с помощью статистического метода. Этот метод основан на том, что свойства макроскопической системы в конечном счете определяются свойствами частиц системы...
74782. Понятие идеального газа. Давление. Основное уравнение молекулярно-кинетической теории газов. Основное уравнение молекулярно-кинетической теории газов 85 KB
  Для вывода основного уравнения молекулярно-кинетической теории рассмотрим одноатомный идеальный газ. Предположим, что молекулы газа движутся хаотически, число взаимных столкновений между молекулами газа пренебрежимо мало по сравнению с числом ударов о стенки сосуда...
74783. Внутренняя энергия системы. Внутренняя энергия идеального газа. Первое начало термодинамики. Примеры 35.5 KB
  Таким образом, можно говорить о двух формах передачи энергии от одних тел к другим: в форме работы и в форме теплоты. Энергия механического движения может превращаться в энергию теплового движения и наоборот.
74784. Теплоемкость (полная, удельная, молярная). Теплоемкость идеального газа (при постоянном давление и объеме). Формула Майера 46.5 KB
  Выражение (53.6) называется уравнением Майера; оно показывает, что Ср всегда больше СV на величину молярной газовой постоянной. Это объясняется тем, что при нагревании газа при постоянном давлении требуется еще дополнительное количество теплоты на совершение работы расширения газа...
74785. Первое начало термодинамики. Круговые, обратимые и необратимые процессы. Тепловая машина Карно и ее кпд 54 KB
  Внутренняя энергия системы может изменяться в результате различных процессов например совершения над системой работы или сообщения ей теплоты. С другой стороны температуру газа и его внутреннюю энергию можно увеличить за счет сообщения ему некоторого количества теплоты...