Показать сообщение отдельно

  #3  
Старый 01.02.2009, 17:59
eLWAux
Постоянный
Регистрация: 15.06.2008
Сообщений: 941
С нами: 9423746

Репутация: 2399


По умолчанию

опять задача из екзамена:
Цитата:
Решить систему уравнений стандартным методом Ньютона без вращения матрицы Якоби(систему линейных уравнений алгебраизма развязать методом Гаусса с формой столбца расписания без выбора главного элемента; якобиан обч. разностным методом), выбирая за начальные приближения
С:
Код:
#include <stdio.h>
#include <math.h>
#include <conio.h>
#define n 2
void main(void)
{
clrscr();
int i, j, L, k, m;
const double E=1e-5;
double x[n], Xprev[n], F[n], J[n][n], a[n][n+1], deltaX[n];
double f(int i,double x[n]);
double Jacobi(int i,int j,double x[n]);
double max(double x[n],double Xprev[n]);

x[0]=0.00000001;     //pochatkove
x[1]=0.00000001;     //nablyzhennya
m=0;

do
{
  for(i=0;i<=n-1;i++)
    Xprev[i]=x[i];
  //----------FORMUYEMO F TA J------------------------
  for(i=0;i<=n-1;i++)
  {
   F[i]=f(i,x);
   for(j=0;j<=n-1;j++)
    J[i][j]=Jacobi(i,j,x);
  }
  //--------------------------------------------------
  for(i=0;i<=n-1;i++)
   for(j=0;j<=n-1;j++)
    a[i][j]=J[i][j];
  a[0][n]=F[0];
  a[1][n]=F[1];

   //------------- GAUSS PO STOVPTSYAH----------------------
  for(L=0;L<=n-2;L++)
  {
    for(k=L+1;k<=n;k++)
    {
      a[L][k]=-a[L][k]/a[L][L];
      for(i=L+1;i<=n-1;i++)
        a[i][k]=a[i][k]+a[i][L]*a[L][k];
    }
  }

    // ------ZVOROTNIY HID-------
  deltaX[n-1]=-a[n-1][n]/a[n-1][n-1];
  for(i=n-2;i>=0;i--)
  {
    deltaX[i]=a[i][n];
    for(k=i+1;k<=n-1;k++)
      deltaX[i]+=a[i][k]*deltaX[k];
  }

  for(i=0;i<=n-1;i++)
    x[i]+=deltaX[i];    //--UTOCHNENE ZNACHENNYA x[i]--
  m++;
}
while(max(x,Xprev)>E);

printf("\t ---ROZVYAZOK SYSTEMY---\n\t|   (ZA %d ITERATSII):   |\n",m);
for(i=0;i<=n-1;i++)
  printf("\t|    x[%d]= %lf    |\n",i,x[i]);
printf("\t -----------------------");
printf("\nPEREVIRKA_1!!!  %lf\n",-x[0]+x[0]*x[0]-x[1]*x[1]-0.1);
printf("PEREVIRKA_2!!!  %lf\n",-x[1]+2*x[0]*x[1]-0.1);

getch();
}

double f(int i,double x[n])
{
 if(i==0)
  return -x[0]+x[0]*x[0]-x[1]*x[1]-0.1;
 else
  return -x[1]+2*x[0]*x[1]-0.1;
}
double Jacobi(int i,int j,double x[n])
{
 const double h=1e-9;
 double prev=f(i,x);
 x[j]+=h;
 return (f(i,x)-prev)/h;
}
double max(double x[n],double Xprev[n])
{
int i;
double z,q;
z=fabs((x[0]-Xprev[0])/Xprev[0]);
for(i=1;i<=n-1;i++)
{
 q=fabs((x[i]-Xprev[i])/Xprev[i]);
 if(q>z)
  z=q;
}
return z;
}
 
Ответить с цитированием