
01.02.2009, 17:59
|
|
Постоянный
Регистрация: 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;
}
|
|
|