数模论坛

 找回密码
 注-册-帐-号
搜索
热搜: 活动 交友 discuz
查看: 3464|回复: 2

一个小程序

[复制链接]
发表于 2003-12-6 22:27:06 | 显示全部楼层 |阅读模式
/*   
                        Graph.h
*/
#include <vector>
#include <complex>
#ifndef _GRAPH_H_
#define _GRAPH_H_

using namespace std;

class Arcc
{
  public:
    Arcc(double r,double c,double l,int i,int j):R(r),C(c),L(l),To1(i),To2(j)
    {}
    Arcc():R(0),C(0),L(0),To1(-1),To2(-1)
    {}
    double R;
    double C;
    double L;
    int To1;
    int To2;
};

class Graph
{
  public:
    Graph(int,bool);
    void AddArcc(Arcc);
    std::complex<double> CaculateImage(double);
    vector<double> CalculateFrequency();
    void SetDelta(double d)
    {
          Delta=d;
    }
    double dF(double w);
    double d2F(double w);
    double d3F(double w);
    bool Judge()
    {
          return JJ;
    }
    vector<complex<double> > Gauss(vector<vector<complex<double> > > A,vector<complex<double> > B);
  protected:
    std::vector<std::vector<int> > Node;//存节点信息的矩阵
    std::vector<Arcc> Arcct;//存边信息的向量
    bool JJ;
    bool output;
    double Delta;
};

Graph::Graph(int i,bool o):Node(i),JJ(false),output(o)
{
      for(int k=0;k<i;++k)
      {
            Node[k].push_back(0);
            Node[k].push_back(k);
      }
}

void Graph::AddArcc(Arcc a)
{
     Arcct.push_back(a);
     ++Node[a.To1][0];
     ++Node[a.To2][0];
}

vector<complex<double> > Graph::Gauss(vector<vector<complex<double> > > A,vector<complex<double> > B)
{
     for(int k=0;k<(int)A.size();++k)
     {
            int g=k;
            for(int i=k;i<(int)A.size();++i)
            {
                  if(abs(A[k])>abs(A[g][k]))
                  {
                        g=i;
                  }
            }
            for(int i=0;i<(int)A.size();++i)
            {
                  complex<double> temp=A[k];
                  A[k]=A[g];
                  A[g]=temp;
            }
            complex<double> temp=B[k];
            B[k]=B[g];
            B[g]=temp;
            if(abs(A[k][k])!=0)
            {
                  complex<double> x=A[k][k];
                  for(int i=0;i<(int)A.size();++i)
                  {
                         A[k]=(A[k]/x);
                  }
                  B[k]=(B[k]/x);
                  for(int i=0;i<(int)A.size();++i)
                  {
                         if(i!=k)
                         {
                               x=A[k];
                               for(int j=0;j<(int)A.size();++j)
                               {
                                       A[j]=A[j]-A[k][j]*x/A[k][k];
                               }
                               B=(B-((B[k]*x)/A[k][k]));
                         }
                  }
            }
     }
     return B;
}

complex<double> Graph::CaculateImage(double w)
{
     vector<vector<int> > node(Node);    //节点转存
     vector<complex<double> > data;  //边权转存
     vector<Arcc> AA(Arcct);  //边转存
     for(int i=0;i<(int)AA.size();++i)
     {
          complex<double> rr(AA.R,0);
          complex<double> ll(0,w*(AA.L));
          if(w*AA.C!=0)
          {
               complex<double> cc(0,AA.C*w);
               data.push_back(rr+ll+1.0/cc);
          }
          else
          {
               data.push_back(rr+ll);
          }
     }
     for(int i=0;i<(int)data.size();++i)
     {
          if(abs(data)==0)
          {
               if(!(((AA.To1==1)&&(AA.To2==0))||((AA.To1==0)&&(AA.To2==1))))
               {
                     int d;
                     vector<complex<double> > :: iterator it1=data.begin();
                     vector<Arcc> :: iterator it2=AA.begin();
                     for(;it1!=data.end();++it1,++it2)
                     {
                            if(((it2->To1==AA.To1)&&(it2->To2==AA.To2))||((it2->To2==AA.To1)&&(it2->To1==AA.To2)))
                            {
                                    it1=data.erase(it1)-1;
                                    it2=AA.erase(it2)-1;
                            }
                     }
                     if(AA.To1>AA.To2)
                     {
                            d=AA.To2;
                            vector<std::vector<int> >:: iterator it3;
                            for(it3=node.begin();it3!=node.end();++it3)
                            {
                                   if((*it3)[1]==AA.To1)
                                   {
                                          it3=node.erase(it3);
                                   }
                                   if((*it3)[1]>AA.To1)
                                   {
                                          --(*it3)[1];
                                   }
                            }
                            for(int k=0;k<(int)data.size();++k)
                            {
                                   if(AA[k].To1==AA.To1)
                                   {
                                          AA[k].To1=d;
                                   }
                                   if(AA[k].To2==AA.To1)
                                   {
                                          AA[k].To2=d;
                                   }
                                   if(AA[k].To1>AA.To1)
                                   {
                                          --AA[k].To1;
                                   }
                                   if(AA[k].To2>AA.To1)
                                   {
                                          --AA[k].To2;
                                   }
                            }
                     }else
                     {
                            d=AA.To1;
                            vector<std::vector<int> >:: iterator it3;
                            for(it3=node.begin();it3!=node.end();++it3)
                            {
                                   if((*it3)[1]==AA.To2)
                                   {
                                         it3=node.erase(it3);
                                   }
                                   if((*it3)[1]>AA.To2)
                                   {
                                         --(*it3)[1];
                                   }
                            }
                            for(int k=0;k<(int)data.size();++k)
                            {
                                   if(AA[k].To1==AA.To2)
                                   {
                                         AA[k].To1=d;
                                   }
                                   if(AA[k].To2==AA.To2)
                                   {
                                         AA[k].To2=d;
                                   }
                                   if(AA[k].To1>AA.To2)
                                   {
                                         --AA[k].To1;
                                   }
                                   if(AA[k].To2>AA.To2)
                                   {
                                         --AA[k].To2;
                                   }
                            }
                     }
               }
          }
     }
     complex<double> zero(0,0);
     complex<double> one(1.0,0);
     vector<vector<complex<double> > > A(node.size()-2);
     vector<complex<double> > B(A.size(),zero);
     for(int i=0;i<(int)A.size();++i)
     for(int j=0;j<(int)A.size();++j)
     {
          A.push_back(zero);
     }
     vector<Arcc> :: iterator it1=AA.begin();
     vector<complex<double> > :: iterator it3=data.begin();
     complex<double> out(0,0);
     int kk(0);
     int kt(0);
     for(int i=0;i<(int)data.size();++i,++it1,++it3)
     {
          if(((it1->To1==1)&&(it1->To2==0))||((it1->To2==1)&&(it1->To1==0)))
          {
               it1=AA.erase(it1)-1;
               if(abs(data)!=0)
               {
                     out=out+(1.0/data);
                     ++kk;
               }
               else
               {
                     ++kt;
               }
               it3=data.erase(it3)-1;

          }
     }
     if(kk!=0)
     {
          out=(1.0/out);
     }
     for(int i=0;i<(int)data.size();++i)
     {
          if((AA.To1==0)||(AA.To2==0))
          {
               if(AA.To1==0)
               {
                     A[AA.To2-2][AA.To2-2]=A[AA.To2-2][AA.To2-2]-(1.0/data);
               }else
               {
                     A[AA.To1-2][AA.To1-2]=A[AA.To1-2][AA.To1-2]-(1.0/data);
               }
          }else if((AA.To1==1)||(AA.To2==1))
          {
               if(AA.To1==1)
               {
                     A[AA.To2-2][AA.To2-2]=A[AA.To2-2][AA.To2-2]-(1.0/data);
                     B[AA.To2-2]=B[AA.To2-2]-(one/data);
               }else
               {
                     A[AA.To1-2][AA.To1-2]=A[AA.To1-2][AA.To1-2]-(1.0/data);
                     B[AA.To1-2]=B[AA.To1-2]-(one/data);
               }
          }else
          {
               A[AA.To2-2][AA.To2-2]=A[AA.To2-2][AA.To2-2]-(1.0/data);
               A[AA.To2-2][AA.To1-2]=A[AA.To2-2][AA.To1-2]+(1.0/data);

               A[AA.To1-2][AA.To1-2]=A[AA.To1-2][AA.To1-2]-(1.0/data);
               A[AA.To1-2][AA.To2-2]=A[AA.To1-2][AA.To2-2]+(1.0/data);
          }
     }
     vector<complex<double> > U(Gauss(A,B));
     complex<double> line(0,0);
     for(int i=0;i<(int)data.size();++i)
     {
          if((AA.To1==1)||(AA.To2==1))
          {
               if(AA.To1==1)
               {
                     line+=((one-U[AA.To2-2])/data);
               }else
               {
                     line+=((one-U[AA.To1-2])/data);
               }
          }
     }
     if(output==false)
     {
          if(kt!=0)
          {
               out=(one/line);
          }else if(kk!=0)
          {
               out=(one/line)+out;
          }
     }else
     {
          if(kk!=0)
          {
               out=(1.0/((line/one)+(1.0)/out));
          }else
          {
               out=(one/line);
          }
     }
     return out;
}

double Graph::dF(double w)
{
     double delta=0.0000000000001;
     complex<double> data=CaculateImage(w);
     complex<double> data0=CaculateImage(w-delta);
     complex<double> data1=CaculateImage(w+delta);
     double x1=(data.imag()-data0.imag())/delta;
     double x2=(data1.imag()-data.imag())/delta;
     return (x1+x2)/2;
}

double Graph::d2F(double w)
{
     double delta=0.0000000000001;
     double data0=dF(w-delta);
     double data=dF(w);
     double data1=dF(w+delta);
     double x0=(data-data0)/delta;
     double x1=(data1-data)/delta;
     return (x0+x1)/2;
}

double Graph::d3F(double w)
{
     double delta=0.0000000000001;
     double data0=d2F(w-delta);
     double data=d2F(w);
     double data1=d2F(w+delta);
     double x0=(data-data0)/delta;
     double x1=(data1-data)/delta;
     return (x0+x1)/2;
}

vector<double> Graph::CalculateFrequency()
{
     JJ=false;
     vector<double> result;
     double Delta1(1000000);
     complex<double> data0=CaculateImage(0);
     double x2(Delta1);
     double x3(0-Delta1);
     complex<double> data2=CaculateImage(x2);
     complex<double> data3=CaculateImage(x3);
     while((data0.imag()*data2.imag()>0)&&(data0.imag()*data3.imag()>0))
     {
          x2+=Delta1;
          x3-=Delta1;
          data2=CaculateImage(x2);
          data3=CaculateImage(x3);
     }
     if(data0.imag()*data3.imag()<=0)
     {
          JJ=true;
     }
     while(data0.imag()*data2.imag()>0)
     {
          x2+=Delta1;
          data2=CaculateImage(x2);
     }
     x3=x2-Delta1;
     while(fabs(x2-x3)>Delta)
     {
          data3=CaculateImage((x2+x3)/2);
          if(data3.imag()*data2.imag()>0)
          {
                x2=(x2+x3)/2;
          }
          else
          {
                x3=(x2+x3)/2;
          }
     }
     result.push_back((x2+x3)/2);
     if(JJ==false)
     {
          x2=x3;
          data2=CaculateImage(x2);
          data3=CaculateImage(x3);
          while(!((data2.imag()*data3.imag()<=0)||(d2F(x3)*(data3.imag()>=0))||(d3F(x3)*(data3.imag()>=0))))
          {
                x3+=Delta1;
          }
          if(data2.imag()*data3.imag()<=0)
          {
                x2=x3-Delta1;
                while(fabs(x2-x3)>Delta)
                {
                       data2=CaculateImage(x3);
                       data3=CaculateImage((x2+x3)/2);
                       if(data3.imag()*data2.imag()>0)
                       {
                               x3=(x2+x3)/2;
                       }
                       else
                       {
                               x2=(x2+x3)/2;
                       }
                }
                result.push_back((x3+x2)/2);
                JJ=true;
          }
     }
     return result;
}

#endif
 楼主| 发表于 2003-12-6 22:27:59 | 显示全部楼层


/*
                     Main.cpp
*/
//---------------------------------------------------------------------------
#include <vcl.h>
#include <complex>
#include <fstream>
#include <process.h>
#include <vector>
#include "Graph.h"
#pragma hdrstop

#include "Main.h"
//---------------------------------------------------------------------------
#pragma resource "*.dfm"
using namespace std;
TMainForm *MainForm;
Graph *g=NULL;
//---------------------------------------------------------------------------
__fastcall TMainForm::TMainForm(TComponent* Owner)
        : TForm(Owner)
{
        Run->Enabled=false;
        Button3->Enabled=true;
        RadioButton2->Checked=true;
}
//----------------------------------------------------------------------------

void __fastcall TMainForm::Button1Click(TObject *Sender)
{
       if((NodeNum->Text!="")&&(g==NULL))
       {
              if(CheckBox1->Checked==true)
              {
                    g=new Graph(NodeNum->Text.ToInt(),true);
              }else
              {
                    g=new Graph(NodeNum->Text.ToInt(),false);
              }
              Button1->Enabled=false;
              CheckBox1->Enabled=false;
              ofstream out("data.txt");
              out<<"电路图中有"<<NodeNum->Text.ToInt()<<"个节点。"<<endl;
              if(CheckBox1->Checked==true)
              {
                    out<<"为二端口电路,节点0和节点1为端口。"<<endl;
              }
              out<<"--------------------------------------------------"<<endl;
              out<<"R(Ω)\tC(F)\tL(H)\tFrom(No.)\tTo(No.)"<<endl;
              AnsiString s1("电路中节点数已经设定好了!请向电路中添加数据!\n");
              AnsiString s2("电路中节点编号从0到");
              AnsiString s3(NodeNum->Text.ToInt()-1);
              MessageDlg(s1+s2+s3,mtInformation,TMsgDlgButtons() << mbOK, 0);
              NodeNum->ReadOnly=true;
              NodeNum->Color=cl3DLight;
       }else if(NodeNum->Text=="")
       {
              MessageDlg("请输入节点中的节点数!",mtInformation,TMsgDlgButtons() << mbOK, 0);
       }
}
//---------------------------------------------------------------------------

void __fastcall TMainForm::AddClick(TObject *Sender)
{
       if((Begin->Text.ToDouble()>=NodeNum->Text.ToDouble())||(End->Text.ToDouble()>=NodeNum->Text.ToDouble())||(Begin->Text.ToDouble()<0)||(End->Text.ToDouble()<0))
       {
              MessageDlg("您所添加的该边的两个端点中有非该图中的节点!",mtError,TMsgDlgButtons() << mbOK, 0);
       }else if(Begin->Text==""||End->Text==""||R->Text==""||L->Text==""||C->Text=="")
       {
              MessageDlg("输入错误!输入数据不可以为空!",mtError,TMsgDlgButtons() << mbOK, 0);
       }else if(Begin->Text.ToDouble()==End->Text.ToDouble())
       {
              MessageDlg("边的起点和终点不可以相同!",mtError,TMsgDlgButtons() << mbOK, 0);
       }
       else if(g)
       {
              Arcc aa(R->Text.ToDouble(),C->Text.ToDouble(),L->Text.ToDouble(),Begin->Text.ToInt(),End->Text.ToInt()) ;
              g->AddArcc(aa);
              ofstream out("data.txt",std::ios:ut|std::ios::app);
              out<<R->Text.ToDouble()<<"\t"<<C->Text.ToDouble()<<"\t"<<L->Text.ToDouble()<<"\t"<<Begin->Text.ToInt()<<"\t\t"<<End->Text.ToInt()<<endl;
              MessageDlg("数据已经添加进去了,请继续添加数据!",mtInformation,TMsgDlgButtons() << mbOK, 0);
              R->Text="";
              C->Text="";
              L->Text="";
              Begin->Text="";
              End->Text="";
       }else if(!g)
       {
              MessageDlg("请先确定图中节点的个数!",mtError,TMsgDlgButtons() << mbOK, 0);
       }
}
//---------------------------------------------------------------------------

void __fastcall TMainForm::Button3Click(TObject *Sender)
{
       if(Result->Text=="")
       {
              MessageDlg("请输入频率!",mtError,TMsgDlgButtons() << mbOK, 0);
       }
       else if(g)
       {
              complex<double> dd(g->CaculateImage(Result->Text.ToDouble()));
              AnsiString s1(dd.real());
              AnsiString s3(dd.imag());
              AnsiString s2("+");
              AnsiString s4("i");
              if(dd.imag()>=0)
              {
                     ww->Text=s1+s2+s3+s4;
              }else if(dd.imag()<0)
              {
                     ww->Text=s1+s3+s4;
              }
              MessageDlg("计算完成!",mtInformation,TMsgDlgButtons() << mbOK, 0);
       }
       else
       {
              MessageDlg("请先输入图中信息!",mtError,TMsgDlgButtons() << mbOK, 0);
       }
}
//---------------------------------------------------------------------------

void __fastcall TMainForm::close(TObject *Sender, TCloseAction &Action)
{
       if(g)
       {
              delete g;
       }
       MessageDlg(" 谢谢使用!\n开发者:倪冉\nKnight-Studio",mtInformation,TMsgDlgButtons() << mbOK, 0);
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
void __fastcall TMainForm::FormCreate(TObject *Sender)
{
       MessageDlg("欢迎使用电路Calculater2.0!\n开发者:倪冉\nKnight-Studio",mtInformation,TMsgDlgButtons() << mbOK, 0);
}
//---------------------------------------------------------------------------

void __fastcall TMainForm::RadioButton1Click(TObject *Sender)
{
       if(RadioButton1->Checked==true)
       {
              Run->Enabled=true;
              Button3->Enabled=false;
              ww->Color=cl3DLight;
              ww->ReadOnly=true;
              delta->ReadOnly=false;
              delta->Color=clWindow;
              Button4->Enabled=true;
              Button5->Enabled=true;
       }
}
//---------------------------------------------------------------------------

void __fastcall TMainForm::RadioButton2Click(TObject *Sender)
{
       if(RadioButton2->Checked==true)
       {
              Run->Enabled=false;
              Button3->Enabled=true;
              ww->ReadOnly=false;
              ww->Color=clWindow;
              delta->ReadOnly=true;
              delta->Color=cl3DLight;
              Button4->Enabled=false;
              Button5->Enabled=false;
       }
}
//---------------------------------------------------------------------------


void __fastcall TMainForm::Button2Click(TObject *Sender)
{
       if(g)
       {
               delete g;
               g=NULL;
               NodeNum->ReadOnly=false;
               NodeNum->Color=clWindow;
               NodeNum->Text="";
               Button1->Enabled=true;
               CheckBox1->Enabled=true;
               Button4->Enabled=true;
               Run->Enabled=false;
               Button3->Enabled=true;
               RadioButton2->Checked=true;
               delta->Text="";
       }

}
//---------------------------------------------------------------------------

void __fastcall TMainForm::RunClick(TObject *Sender)
{
       if(g)
       {
              vector<double> data=g->CalculateFrequency();
              //Result->Text=g->CalculateFrequency();
              if((g->Judge())&&(data.size()==1))
              {
                      Result->Text=data[0];
                      MessageDlg("该解为唯一正解!",mtInformation,TMsgDlgButtons()<<mbOK,0);
              }
              else if((g->Judge())&&(data.size()==2))
              {
                      AnsiString s1(data[0]);
                      AnsiString s3(data[2]);
                      AnsiString s2("OR");
                      Result->Text=s1+s2+s3;
                      MessageDlg("存在两个可以求出的正解!",mtInformation,TMsgDlgButtons()<<mbOK,0);
              }
              else
              {
                      Result->Text=data[0];
                      MessageDlg("该解为一个正解,可能还有其他的正解,不能确定!",mtInformation,TMsgDlgButtons()<<mbOK,0);
              }
       }
       else
       {
              MessageDlg("请先输入电路图中的信息!",mtError,TMsgDlgButtons()<<mbOK,0);
       }
}
//---------------------------------------------------------------------------


void __fastcall TMainForm::Button4Click(TObject *Sender)
{
       if(g)
       {
              g->SetDelta(delta->Text.ToDouble());
              delta->ReadOnly=true;
              delta->Color=cl3DLight;
              Button4->Enabled=false;
       }
       else
       {
              MessageDlg("请先输入电路图中的信息!",mtError,TMsgDlgButtons()<<mbOK,0);
       }
}
//---------------------------------------------------------------------------

void __fastcall TMainForm::Button5Click(TObject *Sender)
{
       delta->ReadOnly=false;
       delta->Color=clWindow;
       delta->Text="";
       Button4->Enabled=true;
}
//---------------------------------------------------------------------------

/*
              main.h
*/
//---------------------------------------------------------------------------
#ifndef MainH
#define MainH
//---------------------------------------------------------------------------
#include <sysutils.hpp>
#include <windows.hpp>
#include <messages.hpp>
#include <sysutils.hpp>
#include <classes.hpp>
#include <graphics.hpp>
#include <controls.hpp>
#include <forms.hpp>
#include <dialogs.hpp>
#include <stdctrls.hpp>
#include <buttons.hpp>
#include <extctrls.hpp>
#include <menus.hpp>
#include <Classes.hpp>
#include <Controls.hpp>
#include <StdCtrls.hpp>
#include <Buttons.hpp>
//---------------------------------------------------------------------------
class TMainForm : public TForm
{
__published:
        TGroupBox *GroupBox1;
        TEdit *NodeNum;
        TLabel *Label1;
        TGroupBox *GroupBox2;
        TLabel *Label2;
        TEdit *Begin;
        TLabel *Label3;
        TEdit *End;
        TLabel *Label4;
        TEdit *R;
        TLabel *Label5;
        TLabel *Label6;
        TButton *Button1;
        TEdit *C;
        TLabel *Label8;
        TEdit *L;
        TLabel *Label9;
        TButton *Add;
        TGroupBox *GroupBox3;
        TLabel *Label10;
        TEdit *Result;
        TButton *Run;
        TLabel *Label7;
        TEdit *ww;
        TButton *Button3;
        TLabel *Label11;
        TGroupBox *GroupBox4;
        TRadioButton *RadioButton1;
        TRadioButton *RadioButton2;
        TLabel *Label12;
        TButton *Button2;
        TGroupBox *GroupBox5;
        TLabel *Label13;
        TEdit *delta;
        TButton *Button4;
        TButton *Button5;
        TCheckBox *CheckBox1;
        TLabel *Label14;
        void __fastcall Button1Click(TObject *Sender);
        void __fastcall AddClick(TObject *Sender);
        void __fastcall Button3Click(TObject *Sender);
        void __fastcall close(TObject *Sender, TCloseAction &Action);
        void __fastcall FormCreate(TObject *Sender);
        void __fastcall RadioButton1Click(TObject *Sender);
        void __fastcall RadioButton2Click(TObject *Sender);
        void __fastcall Button2Click(TObject *Sender);
        void __fastcall RunClick(TObject *Sender);
        void __fastcall Button4Click(TObject *Sender);
        void __fastcall Button5Click(TObject *Sender);

private:        // private user declarations
public:         // public user declarations
        virtual __fastcall TMainForm(TComponent* Owner);
};
//---------------------------------------------------------------------------
extern TMainForm *MainForm;
//---------------------------------------------------------------------------
#endif
/*
                        Calculter2.0.cpp
*/
//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop
//---------------------------------------------------------------------------
USEFORM("Main.cpp", MainForm);
//---------------------------------------------------------------------------
WINAPI WinMain(HINSTANCE, HINSTANCE, LPSTR, int)
{
        try
        {
                 Application->Initialize();
                 Application->CreateForm(__classid(TMainForm), &MainForm);
                 Application->Run();
        }
        catch (Exception &exception)
        {
                 Application->ShowException(&exception);
        }
        catch (...)
        {
                 try
                 {
                         throw Exception("");
                 }
                 catch (Exception &exception)
                 {
                         Application->ShowException(&exception);
                 }
        }
        return 0;
}
//---------------------------------------------------------------------------

 楼主| 发表于 2003-12-6 22:31:13 | 显示全部楼层
以上的代码可在C++Builder6.0中编译
是本人前几天做的一个算电路阻抗和谐振频率的
您需要登录后才可以回帖 登录 | 注-册-帐-号

本版积分规则

小黑屋|手机版|Archiver|数学建模网 ( 湘ICP备11011602号 )

GMT+8, 2024-11-27 04:02 , Processed in 0.054653 second(s), 18 queries .

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表