Форум программистов CodeGuru
18 Январь 2018, 16:13:13 *
Добро пожаловать, Гость. Пожалуйста, войдите или зарегистрируйтесь.

Войти
Новости:
 
   Начало   Помощь Войти Регистрация  
Страниц: [1]   Вниз
  Печать  
Автор Тема: Построение апроксимации эллипсодов  (Прочитано 14234 раз)
0 Пользователей и 1 Гость смотрят эту тему.
gray
Новичок
*
Офлайн Офлайн

Сообщений: 1


Просмотр профиля
« : 28 Июнь 2010, 14:01:54 »

Есть программа для построения 2 эллипсоидов и эллипсоид,  аппроксимирующий  сверху  объединение. Можно ли построить такую эллипсоиду, которая содержала бы 2 эллипсоида?
вот код программы:
Код:
unit Unit1;
{ Дано два эллипсоида:
Ei =(x принадлежит R^n: (x-ai)^t*Qi^(-1)*(x-ai) <=1 ),
с центром в т. ai, i=1,2. Здесь Q1, Q1 > 0 - n x n  симметрические положительно
заданные матрицы.
Требуется:
построить эти эллипсоиды и  их объединение. Построить  также
эллипсоид,  аппроксимирующий  сверху  объединение. }                       
interface
uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, Grids, ExtCtrls, TeeProcs, TeEngine, Chart, Series, opengl;
type
  TForm1 = class(TForm)
    StringGrid1: TStringGrid;
    StringGrid2: TStringGrid;
    Label1: TLabel;
    Edit1: TEdit;
    Button1: TButton;
    StringGrid4: TStringGrid;
    StringGrid5: TStringGrid;
    StaticText1: TStaticText;
    StaticText2: TStaticText;
    StaticText3: TStaticText;
    StaticText4: TStaticText;
    Button6: TButton;
    StringGrid6: TStringGrid;
    StringGrid7: TStringGrid;
    Label2: TLabel;
    Label3: TLabel;
    Label4: TLabel;
    Edit2: TEdit;
    Label5: TLabel;
    Edit3: TEdit;
    Button2: TButton;
    procedure Button1Click(Sender: TObject);
    procedure Button2Click(Sender: TObject);
    procedure Button6Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;
const m=5;
      pi=3.14;
      eps=0.00001; { all numbers less than eps are equal 0 }
type mas=array[1..m,1..m] of real;
type mas2=array[1..m,1..10000] of real;
type vec=array[1..m] of real;
type vec2=array[1..10000] of real;
procedure pixi(n: integer; Q, b: mas;  var A: vec; C1,C2,C3:real);
var n, t, e: integer;
    K1, K2, a, b,c, Q1, Q2, obr: mas;
    A1, A2: vec;
    p1, x1, x2, p2: vec2;
    hF, hG, det, C1,C2,C3: real;
Form1: TForm1;
implementation
uses Unit2;
{$R *.dfm}
procedure Opr(n: integer; Q: mas; var det: real);
//нахождение определителя
var k,j,i: integer;
    r: real; a: mas;
  begin
    for i:=1 to n do
    for j:=1 to n do
    a[i,j]:=Q[i,j];
    det:=1;
    for k:=1 to n do
    begin
      det:=det*a[k,k];
      for j:=k+1 to n do
      begin
        r:=a[j,k]/a[k,k];
        for i:=k to n do
        begin
          a[j,i]:=a[j,i]-r*a[k,i];
        end;
      end;
    end;
  end;
procedure Obratis(n:integer; Q: mas; var b: mas);
var k, i, j: integer;
    c: array[0..30, 0..30] of real;
begin
 for k:=1 to n do
 begin
      for i:=1 to n do
       for j:=1 to n do
       begin
            if (i=k) and (j=k) then
               c[i,j] := 1/b[i,j];
               if (i=k) and (j<>k) then
                  c[i,j] := -b[i,j]/b[k,k];
               if (i<>k) and (j=k) then
                  c[i,j] := b[i,k]/b[k,k];
               if (i<>k) and (j<>k) then
                  c[i,j] := b[i,j] - b[k,j] * b[i,k]/b[k,k];
       end;
      for i:=1 to n do
       for j:=1 to n do b[i, j]:= c[i, j];
 end;
end;
procedure Obratim(n:integer; Q, b: mas; var obr:mas);
//Вычисление обратной матрицы
var i, j: integer;
    i1, j1: integer;
begin
 i1:= 1;
 j1:= 1;
 try
  for i:=1 to n do
    for j:=1 to n do
    begin
      i1:= i;
      j1:= j;
      b[i, j]:= Q[i,j];
    end;
 except
  ShowMessage('Ошибка при вводе числа в столбце ' + IntToStr(j1 + 1) +
  ' строке ' + IntToStr(i1 + 1) + '!');
  exit;
 end;
 Obratis(n,Q,b);
 for i:=1 to n do
  for j:=1 to n do
    obr[i, j]:=Round(1000*b[j, i])/1000;
end;
function Virazh(n: integer; obr: mas; P: vec): real;
//Выполнение вычисления s:=s+p[i]*Q[i,j]^(-1)*p[j], где Q[i,j] - обратная матрица
var s: real;
    i, j: integer;
begin
  s:=0;
  for i:=1 to n do
  for j:=1 to n do
  s:=s+P[i]*obr[i,j]*P[j];
 // if s<0 then s:=s*(-1);
  Virazh:=s;
end;
procedure pixi(n: integer; Q, b: mas;  var A: vec; C1,C2,C3:real);
//нахождение X Y Z и построение эллипсоида
 var
  //adeg,bdeg:integer;
  arad,brad:real;
  P: vec;
begin
 brad:=0;
 arad:=0;
 obratim(n,Q,b,obr);
 glColor3f (C1, C2, C3);
 glPushMatrix;
 glBegin(GL_LINE_STRIP);//Линии вдоль эллипсоида
    //for bdeg:=0 to 180 do
    while brad<6.283185 do
     begin
         //for adeg:=0 to 180 do
         while arad<6.283185 do
             begin
                 //arad:=adeg*2*Pi/180;
                 //brad:=bdeg*2*Pi/180;
                 P[1]:=cos(brad)*cos(arad);
                 P[2]:=cos(brad)*sin(arad);
                 P[3]:=sin(brad);
                 r:=1/sqrt(Virazh(n,obr,P));
                 glVertex3f(A[1]+(P[1]*r),A[2]+(P[2]*r),A[3]+(P[3]*r));
                 arad:=arad+hG;
             end;
     arad:=0;
     brad:=brad+hF;
     end;
 glEnd;
 glBegin(GL_LINE_STRIP);//Линии поперёк эллипсоида
      //for adeg:=0 to 90 do
      while arad<3.141593 do
     begin
         //for bdeg:=0 to 180 do
         while brad<6.283185 do
             begin
                 //arad:=adeg*2*Pi/180;
                 //brad:=bdeg*2*Pi/180;
                 P[1]:=cos(brad)*cos(arad);
                 P[2]:=cos(brad)*sin(arad);
                 P[3]:=sin(brad);
                 r:=1/sqrt(Virazh(n,obr,P));
                 glVertex3f(A[1]+(P[1]*r),A[2]+(P[2]*r),A[3]+(P[3]*r));
                 brad:=brad+hF;
             end;
     arad:=arad+hG;
     brad:=0;
     end;
  glEnd;
  glPopMatrix;
end;
procedure RazmernstrG(n:integer; var strG:TstringGrid; st:string);
//задание размерности таблицы
var i:integer;
  begin
   strG.ColCount:=n+1;
   strG.RowCount:=n+1;
   strG.Cells[0,0]:=st;
   for i:=1 to n do
    begin
     strG.Cells[i,0]:=IntToStr(i);
     strG.Cells[0,i]:=IntToStr(i);
    end;
  end;
procedure RazmernVec(n:integer; var strG:TstringGrid; st:string);
//задание размерности векторов
var i:integer;
  begin
   strG.ColCount:=n+1;
   //strG.RowCount:=n+1;
   strG.Cells[0,0]:=st;
   for i:=1 to n do
    begin
     strG.Cells[i,0]:=IntToStr(i);
     strG.Cells[0,i]:=IntToStr(i);
    end;
  end;
procedure TForm1.Button1Click(Sender: TObject);
//задание размерностей для таблиц и шага для построения эллипсоидов
 begin
  n:=StrToInt(Edit1.Text);
  hF:=StrToFloat(Edit2.Text);
  hG:=StrToFloat(Edit3.Text);
  RazmernStrG(n,stringGrid1,'K1');
  RazmernStrG(n,stringGrid2,'K2');
  RazmernStrG(n,stringGrid4,'Q1');
  RazmernStrG(n,stringGrid5,'Q2');
  RazmernVec(n,stringGrid6,'A1');
  RazmernVec(n,stringGrid7,'A2');
 end;
procedure inpm(n:integer; var K:mas; strG:TstringGrid);
//ввод матрицы а размера n*n из таблицы strG
var i,j:integer;
 begin
  for i:=1 to n do
  for j:=1 to n do
  begin
   //K[i,j]:=StrToInt(strg.Cells[j,i]);
   K[i,j]:=random(10);
   StrG.Cells[j,i]:=FloatToStr(K[i,j]);
  end;
end;
procedure inpv(n:integer; var A:vec; strG:TstringGrid);
//ввод векторов А1, А2
var i: integer;
 begin
   for i:=1 to n do
  //A[i]:=StrToInt(strG.Cells[i,1]);
  begin
   A[i]:=random(3)-1;
   strG.Cells[i,1]:=FloatToStr(A[i]);
  end;
 end;
procedure trnspm(n: integer; var K1,a: mas);
//транспонирование матриц
var i,j: integer;
begin
 for i:=1 to n do
 for j:=1 to n do
 a[i,j]:=K1[j,i];
end;
procedure simmetrm(n: integer; var K1,a,Q1: mas);
//получение положительно определенных симметрических матриц пудем
// умножения A транспонированной на А
var i,j,k: integer; s: real;
begin
 for i:=1 to n do
 for j:=1 to n do
 begin
  s:=0;
  for k:=1 to n do
  s:=s+K1[i,k]*a[k,j];
  Q1[i,j]:=s;
 end;
end;
procedure vivodQ1Q2(n:integer; var Q1:mas; StrG:TStringGrid);
//Вывод положительно определенных симметрических матриц на экран
var i,j:integer;
begin
 simmetrm(n, K1,a,Q1);
 simmetrm(n, K2,b,Q2);
 for i:=1 to n do
 for j:=1 to n do
 StrG.Cells[i,j]:=FloatToStr(Q1[j,i]);
end;
procedure TForm1.Button6Click(Sender: TObject);
//Нахождение Q1, Q2 и их вывод на экран
begin
  inpm(n,K1,stringGrid1);
  inpm(n,K2,stringGrid2);
  inpv(n,A1,StringGrid6);
  inpv(n,A2,StringGrid7);
 trnspm(n, K1,a);
 trnspm(n, K2,b);
 vivodQ1Q2(n,Q1,StringGrid4);
 vivodQ1Q2(n,Q2,StringGrid5);
end;
procedure TForm1.Button2Click(Sender: TObject);
//Переход ко 2ой форме
 begin
  form2.Show;
 end;
end.

и
Код:
unit Unit2;
interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  OpenGL, ExtCtrls,Math, Menus, ComCtrls, StdCtrls,xpman, unit1;
type
  TForm2 = class(TForm)
    procedure FormCreate(Sender: TObject);
    procedure FormPaint(Sender: TObject);
    procedure FormDestroy(Sender: TObject);
    procedure FormResize(Sender: TObject);
    procedure FormMouseDown(Sender: TObject; Button: TMouseButton;
      Shift: TShiftState; X, Y: Integer);
    procedure FormMouseMove(Sender: TObject; Shift: TShiftState; X,
      Y: Integer);
    procedure FormMouseUp(Sender: TObject; Button: TMouseButton;
      Shift: TShiftState; X, Y: Integer);
    procedure FormMouseWheelDown(Sender: TObject; Shift: TShiftState;
      MousePos: TPoint; var Handled: Boolean);
    procedure FormMouseWheelUp(Sender: TObject; Shift: TShiftState;
      MousePos: TPoint; var Handled: Boolean);
  private
    DC : HDC;
    hrc: HGLRC;
    down,myrotate:boolean;
    xm,ym:integer;
    { Private declarations }
  public
    { Public declarations }
  end;
var
  Form2: TForm2;
  x,y,z,arad,brad:real;
  r:real;
Procedure Axes;
implementation
{$R *.dfm}
procedure SetDCPixelFormat (hdc : HDC);
//Задаем формат окна
var
 pfd : TPixelFormatDescriptor;
 nPixelFormat : Integer;
begin
 FillChar (pfd, SizeOf (pfd), 0);
 pfd.dwFlags  := PFD_DRAW_TO_WINDOW or PFD_SUPPORT_OPENGL or PFD_DOUBLEBUFFER;
 nPixelFormat := ChoosePixelFormat (hdc, @pfd);
 SetPixelFormat (hdc, nPixelFormat, @pfd);
end;
procedure TForm2.FormCreate(Sender: TObject);
begin
 dc:=getdc(handle);
 SetDCPixelFormat(DC);
 hrc := wglCreateContext(DC);
 wglMakeCurrent(DC, hrc);
 glloadidentity;
 glClearColor (0.3, 0.0, 10.0,1.0);
 glColor3f (1.5, 9.0, 5.9);
 glEnable (GL_DEPTH_TEST);
 myrotate:=true;
end;
procedure TForm2.FormPaint(Sender: TObject);
//Процедура построения эллипсоида
begin
 glClear (GL_COLOR_BUFFER_BIT or GL_DEPTH_BUFFER_BIT);
   axes;
   pixi(n,Q1,b,A1,0,0,0);
   pixi(n,Q2,b,A2,1,1.32,0.3);
   //pixi(n,Q2,b,A1,0.3,1,0);
 SwapBuffers(DC);
end;
procedure TForm2.FormDestroy(Sender: TObject);
begin
 wglMakeCurrent(0, 0);
 wglDeleteContext(hrc);
 getdc(Handle);
 ReleaseDC (Handle, DC);
 DeleteDC (DC);
end;
procedure TForm2.FormResize(Sender: TObject);
begin
    glViewport(0, 0, ClientWidth, ClientHeight);
    glMatrixMode (GL_PROJECTION);
    glLoadIdentity;
    gluPerspective(30.0, ClientWidth / ClientHeight, 1.0, 15.0);
    glMatrixMode (GL_MODELVIEW);
    glLoadIdentity;
    glTranslatef(0, 0, -10.0);
    glRotatef(30, 1.0, 0.0, 0.0);
    glRotatef(70.0,0.0,1.0,0.0);
    InvalidateRect(Handle, nil, False);
end;
procedure TForm2.FormMouseDown(Sender: TObject; Button: TMouseButton;
  Shift: TShiftState; X, Y: Integer);
begin
    down:=true;
    xm:=x;
    ym:=y
end;
procedure TForm2.FormMouseMove(Sender: TObject; Shift: TShiftState; X,
  Y: Integer);
begin
    if down then begin
    if myrotate then
        begin
            glRotatef (X - Xm, 0.0, 1.0, 0.0);
            glRotatef (Y - Ym, 1.0, 0.0, 0.0);
            InvalidateRect(Handle, nil, False);
            Xm := X;
            Ym := Y;
        end;
    end;
end;
procedure TForm2.FormMouseUp(Sender: TObject; Button: TMouseButton;
  Shift: TShiftState; X, Y: Integer);
begin
    down:=false;
end;
procedure TForm2.FormMouseWheelDown(Sender: TObject; Shift: TShiftState;
  MousePos: TPoint; var Handled: Boolean);
begin
    glscalef(1.2,1.2,1.2);
    InvalidateRect(Handle, nil, False);
end;
procedure TForm2.FormMouseWheelUp(Sender: TObject; Shift: TShiftState;
  MousePos: TPoint; var Handled: Boolean);
begin
    glscalef(1/1.2,1/1.2,1/1.2);
    InvalidateRect(Handle, nil, False);
end;
Procedure Axes;
//Процедура построения Осей координат
begin
glPushMatrix;
glColor3f(1,1,1);
glBegin (GL_LINES);
   glVertex3f(0,0,0);
   glVertex3f(10,0,0);
   glVertex3f(0,0,0);
   glVertex3f(0,10,0);
   glVertex3f(0,0,0);
   glVertex3f(0,0,10);
 glEnd;
glcolor3f(0,0,0);
glBegin (GL_LINES);
glVertex3f (10.1,-0.2, 0.5);
glVertex3f (10.1,0.2, 0.1);
glVertex3f (10.1,-0.2, 0.1);
glVertex3f (10.1,0.2, 0.5);
glEnd;
glBegin (GL_LINES);
glVertex3f (0.0, 10.1, 0.0);
glVertex3f (0.0, 10.1, -0.1);
glVertex3f (0.0, 10.1, 0.0);
glVertex3f (0.1, 10.1, 0.1);
glVertex3f (0.0, 10.1, 0.0);
glVertex3f (-0.1, 10.1, 0.1);
glEnd;
glBegin (GL_LINES);
glVertex3f (0.1, -0.1, 10.1);
glVertex3f (-0.1, -0.1, 10.1);
glVertex3f (0.1, 0.1, 10.1);
glVertex3f (-0.1, 0.1, 10.1);
glVertex3f (-0.1, -0.1, 10.1);
glVertex3f (0.1, 0.1, 10.1);
glEnd;
glPopMatrix;
end;
end.
Записан
3V
Администратор
Ветеран
*****
Офлайн Офлайн

Сообщений: 1347



Просмотр профиля WWW
« Ответ #1 : 30 Июнь 2010, 19:08:31 »

Можно ли построить такую эллипсоиду, которая содержала бы 2 эллипсоида?

Можно, да (инфа 100%!).
Записан

holdmann
Пользователь
***
Офлайн Офлайн

Сообщений: 262



Просмотр профиля
« Ответ #2 : 01 Июль 2010, 00:33:56 »

Можно, да (инфа 100%!).

А где матчасть?
Записан

Елси вы хотите купить, продать, отремонтировать автомобиль в Ижевске: Вам сюда =)
(c)holdmann
Страниц: [1]   Вверх
  Печать  
 
Перейти в:  

Powered by MySQL Powered by PHP Powered by SMF 1.1.21 | SMF © 2015, Simple Machines Valid XHTML 1.0! Valid CSS!