Έστω f(x,y) η πρώτη συνάρτηση και g(x,y) η δεύτερη συνάρτηση και οι μερικές παράγωγοι θf/θx, θf/θy, θg/θx, θg/θy. Στο εξής θα τα συμβολίζω ως f,g,fx,fy,gx,gy.
Στη συνέχεια θα εφαρμόσεις τους αναδρομικούς τύπους:
$x_{n+1}=x_n - \frac{f*gy-g*fy}{fx*gy-gx*fy}$ και $y_{n+1}=y_n - \frac{g*fx-f*gx}{fx*gy-gx*fy}$ με χ0,y0 την αρχική "μαντεψιά" όπως κάνεις και στη Newton για συνάρτηση μιας μεταβλητής.
Το πρόγραμμα θα τερματίζει όταν: f(xn,yn)==g(xn,yn)==abs(eps) (Τότε λύση: (xn,yn) ) ή μετά από τάδε (συνήθως 100) επαναλήψεις (απόκλιση).
Η μέθοδος υπάρχει περίπτωση να βγάλει κάποια overflow και άλλα περίεργα αποτελέσματα, αυτό θα συμβαίνει γιατί για τα starting coords που έδωσες θα έχεις απόκλιση της μεθόδου οπότε δοκίμασε άλλα σημεία. Αν συγκλίνει πάντως θα συγκλίνει γρήγορα (το πολύ πολύ 10 loops).
Ένα παράδειγμα της μεθόδου σε C++- Κώδικας: Επιλογή όλων
#include <iostream.h>
#include <math.h>
void main(){
void Newton2d(double x0, double y0, int max_loops, double tol);
Newton2d( 13,14,100,pow(10,-6) );
}
inline double f(double x,double y) {
double sin(double);
return 1-x-y*sin(3.14159*x/2);
}
inline double f_x(double x,double y){
double cos(double);
return -1-y*cos(3.14159*x/2)*3.14159/2;
}
inline double f_y(double x,double y){
double sin(double);
return -sin(3.14159*x/2);
}
inline double g(double x,double y){
double exp(double),cos(double);
return exp(-y)+10*cos(3.14159*x)-2;
}
inline double g_x(double x,double y){
double sin(double);
return -31.4159*sin(3.14159*x);
}
inline double g_y(double x,double y){
double exp(double);
return -exp(-y);
}
inline double dbl_abs(double num){
return num>0?num:-num;
}
void Newton2d(double x0, double y0, int max_loops, double tol) {
//Variable Declaration
double f(double,double),g(double,double),
f_x(double,double),g_x(double,double),
f_y(double,double),g_y(double,double),
dbl_abs(double),x,y;
int i=1;
//main function
while(i<=max_loops){
x = x0 - ( f(x0,y0)*g_y(x0,y0) - g(x0,y0)*f_y(x0,y0) )/( f_x(x0,y0)*g_y(x0,y0) - g_x(x0,y0)*f_y(x0,y0) );
y = y0 - ( g(x0,y0)*f_x(x0,y0) - f(x0,y0)*g_x(x0,y0) )/( f_x(x0,y0)*g_y(x0,y0) - g_x(x0,y0)*f_y(x0,y0) );
//Test-------------------------------------------------
cout << "Loop "<< i << ":(" << x << "," << y << ").\n";
//-----------------------------------------------------
if( dbl_abs(x-x0) < tol && dbl_abs(y-y0) < tol ){
cout << "The solution is: (" << x <<","<< y <<").\n";
return;
}
++i;
x0=x;
y0=y;
}
cout << "The method has failed after " << max_loops << " loops.\n";
}