Ndtyjky

The following code works like a charm to solve a system of differential equations in it(fcn function in the code), with correct initial values. However, the point of the task is to replace initial values y_1(0) and y_2(0) with some random values, and implement some iterative method to find the correct initial values to solve the equation. I already know how to check if the value is correct value, since by definition output of ddopri 5 should give y_2(1) and y_3(1) as 0. How do I implement Newton Raphson for this problem?

#include<stdio.h>

#include<math.h>

#include<stdbool.h>



double ddopri5(void fcn(double, double *, double *), double *y);

double alpha;

void fcn(double t, double *y, double *f);

double eps;



int main(void){

    double y[4];

    //eps = 1.e-9;

    printf("Enter alpha:n");   

    scanf("%lg", &alpha);

    printf("Enter epsilon:n"); 

    scanf("%lg", &eps);

    y[0]=1.0;//x1(0)

    y[1]=-1.22565282791;//x2(0)

    y[2]=-0.274772807644;//p1(0)

    y[3]=0.0;//p2(0)

    ddopri5(fcn, y);



}





void fcn(double t, double *y, double *f){

/*  double h = 0.25;*/

    f[0] = y[1];

    f[1] = y[3] - sqrt(2)*y[0]*exp(-alpha*t);

    f[2] = sqrt(2)*y[3]*exp(-alpha*t) + y[0];

    f[3] = -y[2];

}





double ddopri5(void fcn(double, double *, double *), double *y){

    double t, h, a, b, tw, chi;

    double w[4], k1[4], k2[4], k3[4], k4[4], k5[4], k6[4], k7[4], err[4], dy[4];

    int i;

    double errabs;

    int iteration;

    iteration = 0;



    //eps = 1.e-9;

    h = 0.1;

    a = 0.0;

    b = 1;//3.1415926535;

    t = a;

    while(t < b -eps){



    printf("%lgn", eps);

        fcn(t, y, k1);

        tw = t+ (1.0/5.0)*h;

        for(i = 0; i < 4; i++){

            /*printf("k1[%i] = %.15lf n", i, k1[i]);*/

            w[i] = y[i] + h*(1.0/5.0)*k1[i];    

        }

        fcn(tw, w, k2);

        tw = t+ (3.0/10.0)*h;

        for(i = 0; i < 4; i++){

            /*printf("k2[%i] = %.15lf n", i, k2[i]);*/

            w[i] = y[i] + h*((3.0/40.0)*k1[i] + (9.0/40.0)*k2[i]);  

        }

        fcn(tw, w, k3);

        tw = t+ (4.0/5.0)*h;

        for(i = 0; i < 4; i++){

            /*printf("k3[%i] = %.15lf n", i, k3[i]);*/

            w[i] = y[i] + h*((44.0/45.0)*k1[i] - (56.0/15.0)*k2[i] + (32.0/9.0)*k3[i]); 

        }

        fcn(tw, w, k4);

        tw = t+ (8.0/9.0)*h;

        for(i = 0; i < 4; i++){

            /*printf("k4[%i] = %.15lf n", i, k4[i]);*/

            w[i] = y[i] + h*((19372.0/6561.0)*k1[i] - (25360.0/2187.0)*k2[i] + (64448.0/6561.0)*k3[i] - (212.0/729.0)*k4[i]);   

        }

        fcn(tw, w, k5);

        tw = t + h;

        for(i = 0; i < 4; i++){

            /*printf("k5[%i] = %.15lf n", i, k5[i]);*/

            w[i] = y[i] + h*((9017.0/3168.0)*k1[i] - (355.0/33.0)*k2[i] + (46732.0/5247.0)*k3[i] + (49.0/176.0)*k4[i] - (5103.0/18656.0)*k5[i]) ;   

        }

        fcn(tw, w, k6);



        tw = t + h;

        for(i = 0; i < 4; i++){

            /*printf("k6[%i] = %.15lf n", i, k6[i]);*/

            w[i] = y[i] + h*((35.0/384.0)*k1[i] + (500.0/1113.0)*k3[i] + (125.0/192.0)*k4[i] - (2187.0/6784.0)*k5[i] + (11.0/84.0)*k6[i]);  

        }

        fcn(tw, w, k7);









        errabs = 0;



        for(i = 0; i < 4; i++){

        /*  printf("k7[%i] = %.15lf n", i, k7[i]);*/

/*          dy[i] = h*((71.0/57600.0)*k1[i]  - (71.0/16695.0)*k3[i] + (71.0/1920.0)*k4[i] - (17253.0/339200.0)*k5[i] + (22.0/525.0)*k6[i]);*/

            dy[i] = h*((35.0/384.0)*k1[i] + (500.0/1113.0)*k3[i] + (125.0/192.0)*k4[i] - (2187.0/6784.0)*k5[i] + (11.0/84.0)*k6[i]);

            /*err[i] =  h*((71.0/57600.0)*k1[i]  + (71.0/16695.0)*k3[i] + (71.0/1920.0)*k4[i] - (17253.0/339200.0)*k5[i] + (22.0/525.0)*k6[i] - (1.0/40.0)*k7[i])*/;

            err[i] =  h*((71.0/57600.0)*k1[i]  - (71.0/16695.0)*k3[i] + (71.0/1920.0)*k4[i] - (17253.0/339200.0)*k5[i] + (22.0/525.0)*k6[i] - (1.0/40.0)*k7[i]);

            /*printf("err[%i] = %.15lf n", i, err[i]);*/

            errabs+=err[i]*err[i];

        }





        errabs = sqrt(errabs);

        printf("errabs = %.15lfn", errabs);

        if( errabs < eps){

            t+= h;

            printf(" FROM IF t t = %.25lf, n h  = %.25lf, n errabs = %.25lf, n iteration = %i . n", t, h, errabs, iteration);

            for(i = 0; i < 4; i++){

                y[i]+=dy[i];            

            }

        } 

/*Avtomaticheskiy vibor shaga*/

        chi=errabs/eps;

        chi = pow(chi, (1.0/6.0));

        if(chi > 10)    chi = 10;

        if(chi < 0.1)   chi = 0.1;

        h*=  0.95/chi;

        if( t + h > b ) h = b - t;      

    /*  for(i = 0; i < 4; i++){

            printf("y[%i] = %.15lf n", i, y[i]);

        }*/



    iteration++;

    printf("t = %.25lf t h = %.25lfn", t, h);

    /*if(iteration > 5) break;*/

    printf("end n");

    for(i = 0; i < 4; i++){

            printf("y[%i] = %.15lf n", i, y[i]);

        }

    if(iteration > 30000) break;

    }

/*  for(i = 0; i < 4; i++){

        printf("y[%i] = %.15lfn", i, y[i]);

    }*/



    return 0;

}

Try this:

Y0=initial_guess

while (true) {

    F=ddopri(Y0);

    Error=F-F_correct

    if (Error small enough)

        break;

    J=jacobian(ddopri, Y0) // this is the matrix dF/dY0

    Y0=Y0-J^(-1)*Error // here you have to solve a linear system

The Jacobian can be obtained using finite differences, i.e. bump up and down the elements of Y one at a time, compute F, take finite differences.

To be clear, element (i,j) of matrix J is dF_i/dY0_j