//====== Глобально заданная размерность системы int n; //====== Граничные условия double U0, UN; //====== Функция вычисления коэффициентов //====== трехдиагональной матрицы double f() { //====== Разовые начальные установки static int raw = -1, k = -1, col = 0; //====== Сдвигаемся по столбцам col++; //====== k считает все элементы //====== Если начинается новая строка if (++k % n == 0) { col = 0; // Обнуляем столбец raw++; // Сдвигаемся по строкам } //====== Выделяем три диагонали return col==raw ? -2. : col == raw-1 || col==raw+1 ? 1. : 0.; } double fu() { //==== Вычисления вектора правых частей по формуле (5) static double dU = (UN-U0)/(n+1), d = U0; return d += dU; } void main() { //======= Размерность задачи и граничные условия n = 4; U0 = 100.; UN = 0.; //======= Размерность valarray (вся матрица) int nn = n*n; //======= Матрица и два вектора valarray a(nn), u(n), v(n); //======= Генерируем их значения generate (&a[0], &a[nn], f); generate (&u[0], &u[n], fu); out("Initial matrix", a); out("Initial vector", u); //======= Умножение матрицы на вектор for (int i=0; i s = a[slice (i*n, n , 1)]; //======= Умножаем на вектор решений //======= Ответ помещаем в вектор v transform(&s[0], &s[n], &u[0], &v[0], multiplies()); //======= Суммируем вектор, получая //======= i-ю компоненту вектора правых частей cout << "\nb[" << i << "] = " << v.sum(); } cout<<"\n\n"; } for (int i=0, id=0; i()); cout << "\nb[" << i << "] = " << v.sum(); } void CChildView::OnPaint() { CPaintDC dc(this); CGraph(m_Points, "Field Distribution", "x[m]","Field").Draw(&dc); } #pragma once #include "Graph.h" class CChildView : public CWnd { //===== Вспомогательные классы будут пользоваться данными friend class CParamDlg; friend class CGraph; private: //===== Контейнер координат точек графика vector m_Points; //===== Вектор источников и свойств среды (см. f и () vector m_f, m_r; //===== Размерность задачи (см. N) int m_n; //===== Параметры double m_k, // Коэффициент k m_L, // Протяженность расчетной области m_g0, // Коэффициенты, задающие ГУ слева m_d0, m_gn, // Коэффициенты, задающие ГУ справа m_dn; CParamDlg *m_pDlg; // Немодальный диалог параметров public: CChildView(); virtual ~CChildView(); virtual BOOL PreCreateWindow(CREATESTRUCT& cs); //===== Изменение размерности задачи void Resize(); //===== Решение системы методом прогонки void Solve(); protected: afx_msg void OnPaint(); DECLARE_MESSAGE_MAP() }; CChildView::CChildView() { m_n = 200; m_k = -0.0005; m_L = 200.; //====== Слева ГУ первого рода Uo=100 m_g0 = 0.; m_d0 = 100.; m_gn = 0.; m_dn = 0.; Resize(); m_pDlg = 0; } CChildView::~CChildView() { m_Points.clear(); m_f.clear(); m_r.clear(); } void CChildView::Resize() { //===== Число узлов равно N+1 (с учетом 0-го узла) int n = m_n + 1; m_Points.resize(n, CDPoint(0.,0.)); m_f.resize(n, 0.); m_r.resize(n, 1.); } void CChildView::Solve() { Resize(); int n = m_n + 1; //======= Коэффициенты разностных уравнений vector a(n), b(n), c(n); //======= Коэффициенты прогонки vector d(n), e(n); double h = m_L / m_n, // Размер шага вдоль оси x hh = h * h; // Квадрат шага //======= Коэффициенты с 0-м индексом не используются a[0] = 0.; b[0] = 0.; c[0] = 0.; //======= Вычисляем координаты x и коэффициенты уравнений m_Points[0].x = 0.; for (int i=1; i < m_n; i++) { m_Points[i].x = i * h; //======= Смотри формулы (4) a[i] = m_r[i-1]/hh; c[i] = m_r[i]/hh; b[i] = - a[i] - c[i] + m_k; } m_Points[m_n].x = m_L; //======= Прямой ход прогонки d[0] = m_g0; // ГУ слева e[0] = m_d0; double den; for (i=1; i < m_n; i++) { //======= Общий знаменатель den = a[i] * d[i-1] + b[i]; d[i] = -c[i] / den; e[i] = (m_f[i] - a[i] * e[i-1]) / den; } //======= Готовимся к обратному ходу den = 1. - m_gn * d[m_n-1]; //======= Случай некорректно поставленной задачи if (den==0.) { MessageBox("ГУ заданы некорректно","Ошибка",MB_OK); return; } //======= Два последних узла используют ГУ справа //======= Смотри формулы (13) m_Points[m_n-1].y = (e[m_n-1] + m_dn * d[m_n-1])/den; m_Points[m_n].y = (m_dn + m_gn* e[m_n-1])/den; //======= Обратный ход прогонки for (i = m_n-2; i >= 0; i--) m_Points[i].y = d[i] * m_Points[i+1].y + e[i]; Invalidate(); } int CChildView::OnCreate(LPCREATESTRUCT lpCreateStruct) { if (CWnd::OnCreate(lpCreateStruct) == -1) return -1; //======= Решаем систему, определенную по умолчанию Solve(); return 0; } #pragma once class CDPoint { public: //======= Две вещественные координаты точки на плоскости double x, y; //======= Стандартный набор конструкторов и операций CDPoint() { x=0.; y=0.; } CDPoint(double xx, double yy) { x=xx; y=yy; } CDPoint& operator=(const CDPoint& pt) { x = pt.x; y = pt.y; return *this; } CDPoint(const CDPoint& pt) { *this = pt; } }; //===== Вспомогательные данные, характеризующие //===== последовательность координат вдоль одной из осей struct TData { //===== Порядок в нормализованном представлении числа int Power; //===== Флаг оси X bool bX; double //======= Экстремумы Min, Max, //======= Множитель (10 в степени Power) Factor, //======= Шаг вдоль оси (мантисса) Step, //======= Реальный шаг dStep, //======= Первая и последняя координаты (мантиссы) Start, End, //======= Первая и последняя координаты dStart, dEnd; }; //===== Класс, реализующий функции плоского графика class CGraph { public: //===== Данные, характеризующие данные вдоль осей TData m_DataX, m_DataY; //===== Контейнер точек графика vector & m_Points; //===== Текущие размеры окна графика CSize m_Size; //===== Экранные координаты центра окна CPoint m_Center; //===== Заголовок и наименования осей CString m_sTitle, m_sX, m_sY; //===== Перо для рисования CPen m_Pen; //===== Два типа шрифтов CFont m_TitleFont, m_Font; //===== Высота буквы (зависит от шрифта) int m_LH, //===== Толщина пера m_Width; //===== Цвет пера COLORREF m_Clr; //======= Методы для управления графиком CGraph(vector& pt, CString sTitle, CString sX, CString sY); virtual ~CGraph(); //===== Заполнение TData для любой из осей void Scale(TData& data); //===== Переход к логическим координатам точек int MapToLogX (double d); int MapToLogY (double d); //===== Изображение в заданном контексте void Draw (CDC *pDC); //===== Изображение одной линии void DrawLine(CDC *pDC); //===== Подготовка цифровой метки на оси CString MakeLabel(bool bX, double& d); }; #include "StdAfx.h" #include "graph.h" //===== Доля окна, занимаемая графиком #define SCALE_X 0.6 #define SCALE_Y 0.6 //===== Внешняя функция нормировки мантисс шагов сетки void gScale (double span, double& step) { //===== Переменная span определяет диапазон изменения //===== значаний одной из координат точек графика //===== Вычисляем порядок числа, описывающего диапазон int power = int(floor(log10(span))); //===== Множитель (zoom factor) double factor = pow(10, power); //===== Мантисса диапазона (теперь 1 < span < 10) span /= factor; //===== Выбираем стандартный шаг сетки if (span<1.99) step=.2; else if (span<2.49) step=.25; else if (span<4.99) step=.5; else if (span<10.) step= 1.; //===== Возвращаем реальный шаг сетки (step*10^power) step *= factor; } void CGraph::Scale (TData& data) { //===== С пустой последовательностью не работаем if (m_Points.empty()) return; //===== Готовимся искать экстремумы data.Max = data.bX ? m_Points[0].x : m_Points[0].y; data.Min = data.Max; //===== Поиск экстремумов for (UINT j=0; j data.Max) data.Max = d; } //===== Максимальная амплитуда двух экстремумов double ext = max(fabs(data.Min),fabs(data.Max)); //===== Искусственно увеличиваем порядок экстремума //===== на 3 единицы, так как мы хотим покрыть 7 порядков, //===== не переходя к экспоненцеальной форме чисел double power = ext > 0.? log10(ext) + 3. : 0.; data.Power = int(floor(power/7.)); //===== Если число не укладывается в этот диапазон if (data.Power != 0) //===== то мы восстанавливаем значение порядка data.Power = int(floor(power)) - 3; //===== Реальный множитель data.Factor = pow(10,data.Power); //===== Диапазон изменения мантиссы double span = (data.Max - data.Min)/data.Factor; //===== Если он нулевой, if (span == 0.) span = 0.5; // то искусственно раздвигаем график //===== Подбираем стандартный шаг для координатной сетки gScale (span, data.Step); //===== Шаг с учетом искусственных преобразований data.dStep = data.Step * data.Factor; //===== Начальная линия сетки должна быть кратна шагу //===== и быть меньше минимума data.dStart = data.dStep * int(floor(data.Min/data.dStep)); data.Start = data.dStart/data.Factor; //===== Вычисляем последнюю линию сетки for (data.End = data.Start; data.End < data.Min/data.Factor + span-1e-10; data.End += data.Step) ; data.dEnd = data.End*data.Factor; } //======= Конструктор класса CGraph CGraph::CGraph (vector& pt, CString sTitle, CString sX, CString sY) : m_Points(pt) { //======= Готовим данные, характеризующие оси координат ZeroMemory(&m_DataX, sizeof(TData)); ZeroMemory(&m_DataY, sizeof(TData)); m_DataX.bX = true; m_DataY.bX = false; m_sTitle = sTitle; m_sX = sX; m_sY = sY; //======= Создаем шрифт для оцифровки осей m_Font.CreateFont(16,0,0,0,100,0,0,0,DEFAULT_CHARSET, OUT_RASTER_PRECIS,CLIP_DEFAULT_PRECIS, DEFAULT_QUALITY,FF_DONTCARE,"Arial"); //======= Выясняем вертикальный размер буквы TEXTMETRIC tm; CClientDC dc(0); dc.SelectObject(&m_Font); dc.GetTextMetrics(&tm); m_LH = tm.tmHeight; //======= Создаем шрифт для вывода заголовка m_TitleFont.CreateFont(24,0,0,0,100,0,0,0,DEFAULT_CHARSET, OUT_RASTER_PRECIS, CLIP_DEFAULT_PRECIS, DEFAULT_QUALITY,FF_DONTCARE,"Times New Roman"); //======= Задаем параметры пера m_Clr = RGB(0,0,255); m_Width = 2; } int CGraph::MapToLogX (double d) { return m_Center.x + int (SCALE_X * m_Size.cx * d); } int CGraph::MapToLogY (double d) { return m_Center.y - int (SCALE_Y * m_Size.cy * d); } //======= Деструктор CGraph::~CGraph(){} void CGraph::Draw(CDC *pDC) { //====== С помощью контекста устройства //====== узнаем адрес окна, его использующего CWnd *pWnd = pDC->GetWindow(); CRect r; pWnd->GetClientRect(&r); //====== Уточняем размеры окна m_Size = r.Size(); m_Center = CPoint(m_Size.cx/2, m_Size.cy/2); int nDC = pDC->SaveDC(); //====== Создаем черное перо для изображения рамки CPen pen(PS_SOLID, 0, COLORREF(0)); pDC->SelectObject(&pen); //====== Преобразуем координаты рамки int lt = MapToLogX(-0.5), rt = MapToLogX(0.5), tp = MapToLogY(0.5), bm = MapToLogY(-0.5); pDC->Rectangle (lt, tp, rt, bm); //====== Задаем цвет и выравнивание текста pDC->SetTextColor(0); pDC->SetTextAlign(TA_LEFT | TA_BASELINE); //====== Выбираем шрифт pDC->SelectObject (&m_Font); //====== Вычисляем атрибуты координатных осей Scale(m_DataX); Scale(m_DataY); //====== Выводим экстремумы функции CString s; s.Format("Min = %.3g",m_DataY.Min); pDC->TextOut(rt+m_LH, tp+m_LH, s); s.Format("Max = %.3g",m_DataY.Max); pDC->TextOut(rt+m_LH, tp+m_LH+m_LH, s); //====== Готовимся изображать координатную сетку CPen gridPen(PS_SOLID, 0, RGB(92,200,178)); pDC->SelectObject(&gridPen); pDC->SetTextAlign(TA_CENTER | TA_BASELINE); //====== Рисуем вертикальные линии сетки for (double x = m_DataX.Start; x < m_DataX.End - m_DataX.Step/2.; x += m_DataX.Step) { //====== Нормируем координату x double xn = (x - m_DataX.Start) / (m_DataX.End - m_DataX.Start) - 0.5; //====== Вычисляем оконную координату int xi = MapToLogX(xn); //====== Пропускаем крайние линии, //====== так как они совпатают с рамкой if (x > m_DataX.Start && x < m_DataX.End) { pDC->MoveTo(xi, bm); pDC->LineTo(xi, tp); } //====== Наносим цифровую метку pDC->TextOut(xi, bm+m_LH, MakeLabel(true, x)); } //===== Повторяем цикл для горизонтальных линий сетки pDC->SetTextAlign(TA_RIGHT | TA_BASELINE); for (double y = m_DataY.Start; y < m_DataY.End - m_DataY.Step/2.; y += m_DataY.Step) { double yn = (y - m_DataY.Start) / (m_DataY.End - m_DataY.Start) - 0.5; int yi = MapToLogY(yn); if (y > m_DataY.Start && y < m_DataY.End) { pDC->MoveTo(lt, yi); pDC->LineTo(rt, yi); pDC->TextOut(lt-m_LH/2, yi, MakeLabel(false, y)); } } //====== Вывод меток осей pDC->TextOut(lt-m_LH/2, tp - m_LH, m_sY); pDC->SetTextAlign(TA_LEFT | TA_BASELINE); pDC->TextOut(rt-m_LH, bm + m_LH, m_sX); //====== Вывод заголовка if (m_sTitle.GetLength() > 40) m_sTitle.Left(40); pDC->SelectObject(&m_TitleFont); pDC->SetTextAlign(TA_CENTER | TA_BASELINE); pDC->TextOut((lt+rt)/2, tp - m_LH, m_sTitle); //====== Вывод линии графика DrawLine(pDC); //====== Восстанавливаем инструменты GDI pDC->RestoreDC(nDC); } void CGraph::DrawLine(CDC *pDC) { //====== Уничтожаем старое перо if (m_Pen.m_hObject) m_Pen.DeleteObject(); //====== Создаем новое m_Pen.CreatePen(PS_SOLID, m_Width, m_Clr); pDC->SelectObject(&m_Pen); double x0 = m_DataX.dStart, y0 = m_DataY.dStart, dx = m_DataX.dEnd - x0, dy = m_DataY.dEnd - y0; for (UINT i=0; i < m_Points.size(); i++) { //====== Нормируем координаты double x = (m_Points[i].x - x0) / dx - .5, y = (m_Points[i].y - y0) / dy - .5; //====== Переход к оконным координатам CPoint pt (MapToLogX(x),MapToLogY(y)); //====== Если точка первая, то ставим перо if (i==0) pDC->MoveTo(pt); else pDC->LineTo(pt); } } CString CGraph::MakeLabel(bool bX, double& v) { CString s = "0.0"; if (v == 0.) return s; //====== Сначала делаем грубую прикидку //====== Пробуем поместиться в 20 позиций s.Format("%20.10f",v); //====== Выясняем порядок шага сетки, //====== переворачивая его знак (трюк) int nDigits = int(ceil(-log10(bX ? m_DataX.Step : m_DataY.Step))); //====== Если все изменения происходят до запятой, //====== то цифры после запятой нас не интересуют if (nDigits <= 0) nDigits = -1; else if (bX) nDigits++; // Эстетическая добавка //====== Слева усекаем все s.TrimLeft(); //====== Справа оставляем минимум разрядов s = s.Left(s.Find(".") + nDigits + 1); int iPower = bX ? m_DataX.Power : m_DataY.Power; //====== Нужен ли порядок? if (iPower != 0) { //====== Нужен, если не поместились в (10^-3, 10^+4) CString add; add.Format("·e%+d",iPower); s += add; } return s; } Элемент Идентификатор Диалог IDD_PARAM Окно редактирования Source: IDC_SOURCE Окно редактирования Start группы Source: IDC_SOURCE1 Окно редактирования End группы Source: IDC_SOURCE2 Окно редактирования Value: IDC_PROP Окно редактирования Start группы Properties: IDC_PROP1 Окно редактирования End группы Properties: IDC_PROP2 Окно редактирования Nodes: IDC_NODES Окно редактирования Distance: IDC_DIST Окно редактирования Decrement: IDC_DECR Окно редактирования g группы Left Boundary: IDC_LEFTG Окно редактирования d группы Left Boundary: IDC_LEFTD Окно редактирования g группы Right Boundary: IDC_RIGHTG Окно редактирования d группы Right Boundary: IDC_RIGHTD Кнопка Add группы Source: IDC_ADDSOURCE Кнопка Add группы Properties: IDC_ADDPROP Кнопка Apply: IDC_APPLY Кнопка Close IDCANCEL #pragma once class CParamDlg : public CDialog { //===== Будем общаться с окном friend class CChildView; DECLARE_DYNAMIC(CParamDlg) public: //===== Будем помнить его адрес CChildView *m_pView; //===== В конструкторе запросим его адрес CParamDlg(CChildView* p); virtual ~CParamDlg(); // Dialog Data enum { IDD = IDD_PARAM }; protected: virtual void DoDataExchange(CDataExchange* pDX); DECLARE_MESSAGE_MAP() }; //==== Интенсивность источника поля double m_Source; //==== Индекс ячейки сетки, где расположено начало источника int m_SrcId1; //==== Индекс ячейки сетки, где расположен конец источника int m_SrcId2; //==== Значение физического свойства ячейки сетки double m_Prop; //==== Индексы начала и конца области со свойством m_Prop int m_PropId1; int m_PropId2; void CParamDlg::DoDataExchange(CDataExchange* pDX) { DDX_Text(pDX, IDC_PROP2, m_PropId2); DDX_Text(pDX, IDC_PROP1, m_PropId1); DDX_Text(pDX, IDC_PROP, m_Prop); DDX_Text(pDX, IDC_SOURCE2, m_SrcId2); DDX_Text(pDX, IDC_SOURCE1, m_SrcId1); DDX_Text(pDX, IDC_SOURCE, m_Source); //========= Обмен с переменными оконного класса DDX_Text(pDX, IDC_NODES, m_pView->m_n); DDX_Text(pDX, IDC_DIST, m_pView->m_L); DDX_Text(pDX, IDC_DECR, m_pView->m_k); DDX_Text(pDX, IDC_LEFTG, m_pView->m_g0); DDX_Text(pDX, IDC_LEFTD, m_pView->m_d0); DDX_Text(pDX, IDC_RIGHTG, m_pView->m_gn); DDX_Text(pDX, IDC_RIGHTD, m_pView->m_dn); CDialog::DoDataExchange(pDX); } void CParamDlg::OnClickedApply(void) { //====== Считываем данные из окон UpdateData(); //====== Заново решаем систему и выводим график m_pView->Solve(); } void CParamDlg::OnClickedAddsource(void) { UpdateData(); //====== Изменяем контейнер m_f (источников поля) for (int i=m_SrcId1; i <= m_SrcId2; i++) { if (0 <= i && i < m_pView->m_n) m_pView->m_f[i] = -m_Source; } m_pView->Solve(); } void CParamDlg::OnClickedAddprop(void) { UpdateData(); //====== Изменяем контейнер m_r (свойств среды) for (int i=m_PropId1; i <= m_PropId2; i++) { if (0 <= i && i < m_pView->m_n && m_Prop > 0.) m_pView->m_r[i] = m_Prop; } m_pView->Solve(); } void CParamDlg::OnClickedCancel(void) { //====== Закрываем немодальный диалог m_pView->m_pDlg = 0; DestroyWindow(); } #include "stdafx.h" #include "Heat.h" #include "ParamDlg.h" IMPLEMENT_DYNAMIC(CParamDlg, CDialog) CParamDlg::CParamDlg(CChildView* p) : CDialog(CParamDlg::IDD, p) { m_pView = p; //===== Начальное значение свойств среды //===== не должно равняться нулю m_Prop = 1.0; m_PropId1 = 0; m_PropId2 = 0; m_Source = 0.0; m_SrcId1 = 0; m_SrcId2 = 0; } CParamDlg::~CParamDlg() { } BOOL CParamDlg::OnInitDialog(void) { CDialog::OnInitDialog(); CRect r; //===== С помощью контекста устройства //===== узнаем размеры всего экрана CClientDC dc(this); int w = dc.GetDeviceCaps(HORZRES); int h = dc.GetDeviceCaps(VERTRES); //===== Узнаем размеры окна диалога GetWindowRect(&r); //===== Смещаем его вправо и вниз r.OffsetRect(w-r.right-10,h-r.bottom-30); MoveWindow(&r); return TRUE; } void CChildView::OnEditParameters(void) { //===== Если диалог не открыт, if (!m_pDlg) { //===== Динамически создаем объект диалогового класса m_pDlg = new CParamDlg(this); //===== и после этого создаем окно диалога m_pDlg->Create(IDD_PARAM); } } void CParamDlg::PostNcDestroy(void) { delete this; }