Расчет двойного интеграла при помощи метода Симпсона
/*****************************************************************************
* .FILE : numeric.c
* .TITLE : Расчет двойного интеграла при помощи метода Симпсона (парабол)
* .DESCR : Курсовой проект по численным методам (1994 год)
* :
* : Выполнил: Петренко В.С. (гр. ПС-301)
* : Проверил: Панюков А.В.
* :
* .NOTE : NOT FOR RENTAL OR SALE.
* : FEDERAL LAW PROVIDES SEVERE CIVIL & CRIMINAL PENALTIES FOR
* : UNAUTHORIZED DUPLICATION OR DISTRIBUTION.
* :
* : (C) '94 by P$P
*****************************************************************************/
#include <math.h>
#include <stdio.h>
/*****************************************************************************
* .NAME : m_Simpson
* .TITLE : Расчет интеграла методом Симпсона (парабол)
* .DESCR :
* :
* .PARAMS : double m_Simpson (double (*func) (double, double),
* : double t_fix, double t_limit, int N);
* : double (*func) (double, double) - подынтегральная ф-я
* : double t_fix - фиксированный первый аргумент
* : double t_limit - верхний предел интегрирования,
* : нижний равен -t_limit
* : int N - число точек разбиения
* .RETURN : Значение вычисленного интеграла
*****************************************************************************/
double m_Simpson (double (*func) (double, double),
double t_fix, double t_limit, int N)
{
double sum1 = 0; /* -¬ */
double sum2 = 0; /* ¦ временные переменные */
double sum3 = 0; /* -- */
double sum; /* конечный результат */
double h = (2 * t_limit) / N; /* шаг сетки */
int i; /* временная */
sum1 = (*func) (t_fix, -t_limit) + (*func) (t_fix, +t_limit);
for (i = 1; i <= N - 1; i++) sum2 += (*func) (t_fix, -t_limit + (i * h));
sum2 *= 2;
for (i = 1; i <= N; i++) sum3 += (*func) (t_fix, -t_limit + ((i-0.5) * h));
sum3 *= 4;
sum = sum1 + sum2 + sum3;
sum = (h / 6) * sum;
return sum;
}
/* */
/* Глобальные переменные */
/* ~~~~~~~~~~~~~~~~~~~~~ */
#define PI 3.1415926536 /* число П */
double k; /* параметр функции - задается пользователем */
int N_MAX; /* число узлов сетки разбиения */
double (*currFunc) (double, double); /* выбранная пользователем функция */
double f1 (double x, double y);
double f2 (double x, double y);
double f3 (double x, double y);
/*****************************************************************************
* .NAME : double F (double dummy, double t)
* .TITLE : Вычисляет внутренний интеграл (G (t)).
* : См. текст курсового проекта.
* .DESCR : П/2
* : -
* : первая вычисляемая функция ¦ G (t + П/2) * sin (t + П/2) dt
* : t -
* : - -П/2
* : где G (t) = ¦ currFunc (t, tau) dtau
* : -
* : -t
* .PARAMS : double F (double dummy, double t);
* : double dummy - фиктивный первый аргумент, при вызове этой
* : функции он не используется, т.к. она
* : функция одного аргумента
* : double - действительный второй аргумент
* .RETURN : Значение функции: G (t) * sin (k * t);
*****************************************************************************/
double F (double dummy, double t)
{
double G;
t = t + PI / 2; /* сдвижка начала координат, чтобы пределы */
/* были симметричны (в нашем случае - на П/2) */
G = m_Simpson (currFunc, t, t, N_MAX);
return G * sin (k * t);
}
/*****************************************************************************
* .NAME : main
* .TITLE : Основная диалоговая функция.
* .DESCR : Запрашивается интересующая пользователя функция,
* : параметр k и число узлов сетки N_MAX.
* : Выводит на экран вычисленное значение интеграла и
* : два справочных значения - П и П/2.
* :
* .PARAMS : void main (void);
* .RETURN :
*****************************************************************************/
void main (void)
{
double integral; /* значение вычисленного интеграла */
int selection; /* номер выбранной функции */
/* массив доступных функций */
double (*functions []) (double, double) = { f1, f2, f3 };
printf ("n Вычисление интеграла методом Симпсона (парабол) ");
printf ("n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ");
printf ("n -- ");
printf ("n I = ¦¦ sin k(x + y) f (x, y) dx dy ");
printf ("n -- ");
printf ("n D ");
printf ("n где D = { (x, y): x, y >= 0; x + y <= П }, f Е C (D)");
printf ("n");
printf ("nДля какой функции рассчитывать: ");
printf ("n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ");
printf ("n 1) f (x, y) = 0.5 * cos (y) ");
printf ("n - -- 0; k != 1 ");
printf ("n ¦ sin x * sin (kx) dx => ¦ ");
printf ("n - L- П/2; k = 1 ");
printf ("n ");
printf ("n 2) f (x, y) = 0.5 - sin (y) ");
printf ("n - ");
printf ("n ¦ x * sin (kx) dx =====> П; k = 1 ");
printf ("n - ");
printf ("n ");
printf ("n 3) f (x, y) = sqrt (x * x + y * y)");
printf ("n");
do
{
printf ("Ваш выбор: ");
scanf ("%d", &selection);
} while (!(1 <= selection && selection <= 3));
printf ("Параметр k: ");
scanf ("%lg", &k);
do
{
printf ("Число узлов сетки N: ");
scanf ("%d", &N_MAX);
} while (!(N_MAX > 0));
printf ("n");
printf ("n Расчет интеграла ...");
currFunc = functions [selection - 1]; /* текущая функция */
integral = m_Simpson (F, 0, PI / 2, N_MAX); /* вычисляем интеграл */
printf ("n Значение интеграла равно: %.12lg", integral); /* вывод */
printf ("n Величины: П = %.12lg; П/2 = %.12lg", PI, PI / 2);
}
/*****************************************************************************
* .FILE : func.c
* .TITLE : Содержит функции пользователя, которые можно
* : изменять без перекомпиляции основной программы
* .DESCR : После изменения этого модуля его необходимо перекомпилировать
* : и слинковать с numeric.obj
* :
* :
* .NOTE : NOT FOR RENTAL OR SALE.
* : FEDERAL LAW PROVIDES SEVERE CIVIL & CRIMINAL PENALTIES FOR
* : UNAUTHORIZED DUPLICATION OR DISTRIBUTION.
* :
* : (C) '94 by P$P
*****************************************************************************/
#include <math.h>
/* выбираемая пользователем функция No.1 */
double f1 (double x, double y)
{ return 0.5 * cos (y); }
/* выбираемая пользователем функция No.2 */
double f2 (double x, double y)
{ return 0.5 - sin (y); }
/* выбираемая пользователем функция No.3 */
double f3 (double x, double y)
{ return sqrt (x * x + y * y); }