// 3元
#include <stdio.h>
double f1(double x, double y, double z)
{
return -y*y * z*z - y*y + 24*y*z - z*z -13;
}
double f2(double x, double y, double z)
{
return -x*x * z*z - x*x + 24*x*z - z*z -13;
}
double f3(double x, double y, double z)
{
return -x*x * y*y - x*x + 24*x*y - y*y -13;
}
double J11(double x, double y, double z) // df1(x,y,z)/dx
{
return 0;
}
double J12(double x, double y, double z) // df1(x,y,z)/dz
{
return -2*z*z*y - 2*y + 24*z;
}
double J13(double x, double y, double z) // df1(x,y,z)/dz
{
return -2*y*y*z + 24*y - 2*z;
}
double J21(double x, double y, double z) // df2(x,y,z)/dx
{
return -2*z*z*x - 2*x + 24*z;
}
double J22(double x, double y, double z) // df2(x,y,z)/dy
{
return 0;
}
double J23(double x, double y, double z) // df2(x,y,z)/dz
{
return -2*x*x*z + 24*x - 2*z;
}
double J31(double x, double y, double z) // df3(x,y,z)/dx
{
return -2*x*y*y - 2*x + 24*y;
}
double J32(double x, double y, double z) // df3(x,y,z)/dy
{
return -2*x*x*y + 24*x - 2*y;
}
double J33(double x, double y, double z) // df3(x,y,z)/dz
{
return 0;
}
void ucalc(double x, double y, double z, double *ux, double *uy, double *uz)
{
double det;
det = J11(x,y,z) * J22(x,y,z) * J33(x,y,z)
+ J21(x,y,z) * J32(x,y,z) * J13(x,y,z)
+ J31(x,y,z) * J12(x,y,z) * J13(x,y,z)
- J11(x,y,z) * J32(x,y,z) * J23(x,y,z)
- J31(x,y,z) * J22(x,y,z) * J13(x,y,z)
- J21(x,y,z) * J12(x,y,z) * J33(x,y,z);
*ux = 1/det * (
(J22(x,y,z)*J33(x,y,z)-J23(x,y,z)*J32(x,y,z))*f1(x,y,z) +
(J13(x,y,z)*J32(x,y,z)-J12(x,y,z)*J33(x,y,z))*f2(x,y,z) +
(J12(x,y,z)*J23(x,y,z)-J13(x,y,z)*J22(x,y,z))*f3(x,y,z));
*uy = 1/det * (
(J23(x,y,z)*J31(x,y,z)-J21(x,y,z)*J33(x,y,z))*f1(x,y,z) +
(J11(x,y,z)*J33(x,y,z)-J13(x,y,z)*J31(x,y,z))*f2(x,y,z) +
(J13(x,y,z)*J21(x,y,z)-J11(x,y,z)*J23(x,y,z))*f3(x,y,z));
*uz = 1/det * (
(J21(x,y,z)*J32(x,y,z)-J22(x,y,z)*J31(x,y,z))*f1(x,y,z) +
(J12(x,y,z)*J31(x,y,z)-J11(x,y,z)*J32(x,y,z))*f2(x,y,z) +
(J11(x,y,z)*J22(x,y,z)-J12(x,y,z)*J21(x,y,z))*f3(x,y,z));
}
int main()
{
double x, y, z, ux, uy, uz;
int i;
x = 1;
y = 1;
z = 1;
for (i = 0; i < 10; i++) {
ucalc(x, y, z, &ux, &uy, &uz);
x -= ux;
y -= uy;
z -= uz;
printf("%d %.20f %.20f %.20f\n", i
, x
, y
, z
); }
return 0;
}
Ly8gM+WFgwojaW5jbHVkZSA8c3RkaW8uaD4KCmRvdWJsZSBmMShkb3VibGUgeCwgZG91YmxlIHksIGRvdWJsZSB6KQp7CglyZXR1cm4gLXkqeSAqIHoqeiAtIHkqeSArIDI0KnkqeiAtIHoqeiAtMTM7Cn0KZG91YmxlIGYyKGRvdWJsZSB4LCBkb3VibGUgeSwgZG91YmxlIHopCnsKCXJldHVybiAteCp4ICogeip6IC0geCp4ICsgMjQqeCp6IC0geip6IC0xMzsKfQpkb3VibGUgZjMoZG91YmxlIHgsIGRvdWJsZSB5LCBkb3VibGUgeikKewoJcmV0dXJuIC14KnggKiB5KnkgLSB4KnggKyAyNCp4KnkgLSB5KnkgLTEzOwp9CmRvdWJsZSBKMTEoZG91YmxlIHgsIGRvdWJsZSB5LCBkb3VibGUgeikJLy8gZGYxKHgseSx6KS9keAp7CglyZXR1cm4gMDsKfQpkb3VibGUgSjEyKGRvdWJsZSB4LCBkb3VibGUgeSwgZG91YmxlIHopCS8vIGRmMSh4LHkseikvZHoKewoJcmV0dXJuIC0yKnoqeip5IC0gMip5ICsgMjQqejsKfQpkb3VibGUgSjEzKGRvdWJsZSB4LCBkb3VibGUgeSwgZG91YmxlIHopCS8vIGRmMSh4LHkseikvZHoKewoJcmV0dXJuIC0yKnkqeSp6ICsgMjQqeSAtIDIqejsKfQpkb3VibGUgSjIxKGRvdWJsZSB4LCBkb3VibGUgeSwgZG91YmxlIHopCS8vIGRmMih4LHkseikvZHgKewoJcmV0dXJuIC0yKnoqeip4IC0gMip4ICsgMjQqejsKfQpkb3VibGUgSjIyKGRvdWJsZSB4LCBkb3VibGUgeSwgZG91YmxlIHopCS8vIGRmMih4LHkseikvZHkKewoJcmV0dXJuIDA7Cn0KZG91YmxlIEoyMyhkb3VibGUgeCwgZG91YmxlIHksIGRvdWJsZSB6KQkvLyBkZjIoeCx5LHopL2R6CnsKCXJldHVybiAtMip4KngqeiArIDI0KnggLSAyKno7Cn0KZG91YmxlIEozMShkb3VibGUgeCwgZG91YmxlIHksIGRvdWJsZSB6KQkvLyBkZjMoeCx5LHopL2R4CnsKCXJldHVybiAtMip4KnkqeSAtIDIqeCArIDI0Knk7Cn0KZG91YmxlIEozMihkb3VibGUgeCwgZG91YmxlIHksIGRvdWJsZSB6KQkvLyBkZjMoeCx5LHopL2R5CnsKCXJldHVybiAtMip4KngqeSArIDI0KnggLSAyKnk7Cn0KZG91YmxlIEozMyhkb3VibGUgeCwgZG91YmxlIHksIGRvdWJsZSB6KQkvLyBkZjMoeCx5LHopL2R6CnsKCXJldHVybiAwOwp9Cgp2b2lkIHVjYWxjKGRvdWJsZSB4LCBkb3VibGUgeSwgZG91YmxlIHosIGRvdWJsZSAqdXgsIGRvdWJsZSAqdXksIGRvdWJsZSAqdXopCnsKCWRvdWJsZQlkZXQ7CgoJZGV0CT0gSjExKHgseSx6KSAqIEoyMih4LHkseikgKiBKMzMoeCx5LHopCgkJKyBKMjEoeCx5LHopICogSjMyKHgseSx6KSAqIEoxMyh4LHkseikKCQkrIEozMSh4LHkseikgKiBKMTIoeCx5LHopICogSjEzKHgseSx6KQoJCS0gSjExKHgseSx6KSAqIEozMih4LHkseikgKiBKMjMoeCx5LHopCgkJLSBKMzEoeCx5LHopICogSjIyKHgseSx6KSAqIEoxMyh4LHkseikKCQktIEoyMSh4LHkseikgKiBKMTIoeCx5LHopICogSjMzKHgseSx6KTsKCSp1eCA9IDEvZGV0ICogKAoJCShKMjIoeCx5LHopKkozMyh4LHkseiktSjIzKHgseSx6KSpKMzIoeCx5LHopKSpmMSh4LHkseikgKwoJCShKMTMoeCx5LHopKkozMih4LHkseiktSjEyKHgseSx6KSpKMzMoeCx5LHopKSpmMih4LHkseikgKwoJCShKMTIoeCx5LHopKkoyMyh4LHkseiktSjEzKHgseSx6KSpKMjIoeCx5LHopKSpmMyh4LHkseikpOwoJKnV5ID0gMS9kZXQgKiAoCgkJKEoyMyh4LHkseikqSjMxKHgseSx6KS1KMjEoeCx5LHopKkozMyh4LHkseikpKmYxKHgseSx6KSArCgkJKEoxMSh4LHkseikqSjMzKHgseSx6KS1KMTMoeCx5LHopKkozMSh4LHkseikpKmYyKHgseSx6KSArCgkJKEoxMyh4LHkseikqSjIxKHgseSx6KS1KMTEoeCx5LHopKkoyMyh4LHkseikpKmYzKHgseSx6KSk7CgkqdXogPSAxL2RldCAqICgKCQkoSjIxKHgseSx6KSpKMzIoeCx5LHopLUoyMih4LHkseikqSjMxKHgseSx6KSkqZjEoeCx5LHopICsKCQkoSjEyKHgseSx6KSpKMzEoeCx5LHopLUoxMSh4LHkseikqSjMyKHgseSx6KSkqZjIoeCx5LHopICsKCQkoSjExKHgseSx6KSpKMjIoeCx5LHopLUoxMih4LHkseikqSjIxKHgseSx6KSkqZjMoeCx5LHopKTsKfQoKaW50IG1haW4oKQp7Cglkb3VibGUJeCwgeSwgeiwgdXgsIHV5LCB1ejsKCWludAlpOwoKCXggPSAxOwoJeSA9IDE7Cgl6ID0gMTsKCWZvciAoaSA9IDA7IGkgPCAxMDsgaSsrKSB7CgkJdWNhbGMoeCwgeSwgeiwgJnV4LCAmdXksICZ1eik7CgkJeCAtPSB1eDsKCQl5IC09IHV5OwoJCXogLT0gdXo7CgkJcHJpbnRmKCIlZCAlLjIwZiAlLjIwZiAlLjIwZlxuIiwgaSwgeCwgeSwgeik7Cgl9CglyZXR1cm4gMDsKfQo=