//------------------- // Solution of Boundary Value Problem with the shooting method. // % Equation to solve: // x''''-(1+t^2)(x'')^2+5x^2=0 // the interval [0,1] // Boundary conditions are: x(0)=1; x'(0)=0; x''(1)=-2; x'''(1)=-3. //-------------------- /* compile with gcc -o NLinearShooting1 NLinearShooting1.c -L/usr/local/dislin -ldislnc -lm -lX11 -lgsl -lgslcblas */ #include "/usr/local/dislin/dislin.h" #include #include #include #include #include #define max 500 int max_steps = 1.; // t in [0, max_steps] int i, j, k, nn; /* Loop counter */ int Kmax=100; float T[max]; float X1[max]; float X2[max]; float X3[max]; float X4[max]; double h = 0.025; // time step double ds=0.01; // step for the derivative // initial conditions for the corresponding i.v.p. double x1 = 1.0; // boundary condition x(0)=1 on the left side of the interval [0, max_steps] double x2 = 0.0; // boundary condition x'(0)=0 on the left side of the interval [0, max_steps] double p = 0.0; // slope to be varied double q = 0.0; // slope to be varied double b3= -2.0; // boundary condition on the right side of the interval [0, max_steps] double b4=-3.0; // boundary condition on the right side of the interval [0, max_steps] double k1x1, k2x1, k3x1, k4x1, k1x2, k2x2, k3x2, k4x2, k1x3, k2x3, k3x3, k4x3, k1x4, k2x4, k3x4, k4x4; /* Function values for RK4 */ double t=0.0; // r.h.s. for the o.d.e's of the 1st order double funct_x1(double t, double x1, double x2, double x3, double x4) { return x2; } double funct_x2(double t, double x1, double x2, double x3, double x4) { return x3; } double funct_x3(double t, double x1, double x2, double x3, double x4) { return x4; } double funct_x4(double t, double x1, double x2, double x3, double x4) { return (1.+pow(t,2))*pow(x3,2)-5.*pow(x1,2); } // RK4 method void rk4(double h, double max_steps, double t, double x1, double x2, double x3, double x4) { T[0] = t; X1[0] = x1; X2[0] = x2; X3[0] = x3; X4[0] = x4; j = 0; while(T[j] < max_steps) { if( (T[j] + h) > max_steps ) h = max_steps - T[j]; /* Calculation of the last step */ t = T[j]; x1 = X1[j]; x2 = X2[j]; x3 = X3[j]; x4 = X4[j]; k1x1 = h * funct_x1(t, x1, x2, x3, x4); k1x2 = h * funct_x2(t, x1, x2, x3, x4); k1x3 = h * funct_x3(t, x1, x2, x3, x4); k1x4 = h * funct_x4(t, x1, x2, x3, x4); k2x1 = h * funct_x1(t+0.5*h, x1 + 0.5 * k1x1, x2 + 0.5 * k1x2, x3 + 0.5 * k1x3, x4 + 0.5 * k1x4); k2x2 = h * funct_x2(t+0.5*h, x1 + 0.5 * k1x1, x2 + 0.5 * k1x2, x3 + 0.5 * k1x3, x4 + 0.5 * k1x4); k2x3 = h * funct_x3(t+0.5*h, x1 + 0.5 * k1x1, x2 + 0.5 * k1x2, x3 + 0.5 * k1x3, x4 + 0.5 * k1x4); k2x4 = h * funct_x4(t+0.5*h, x1 + 0.5 * k1x1, x2 + 0.5 * k1x2, x3 + 0.5 * k1x3, x4 + 0.5 * k1x4); k3x1 = h * funct_x1(t+0.5*h, x1 + 0.5 * k2x1, x2 + 0.5 * k2x2, x3 + 0.5 * k2x3, x4 + 0.5 * k2x4 ); k3x2 = h * funct_x2(t+0.5*h, x1 + 0.5 * k2x1, x2 + 0.5 * k2x2, x3 + 0.5 * k2x3, x4 + 0.5 * k2x4 ); k3x3 = h * funct_x3(t+0.5*h, x1 + 0.5 * k2x1, x2 + 0.5 * k2x2, x3 + 0.5 * k2x3, x4 + 0.5 * k2x4 ); k3x4 = h * funct_x4(t+0.5*h, x1 + 0.5 * k2x1, x2 + 0.5 * k2x2, x3 + 0.5 * k2x3, x4 + 0.5 * k2x4 ); k4x1 = h * funct_x1(t+h, x1 + k3x1, x2 + k3x2, x3 + k3x3, x4 + k3x4); k4x2 = h * funct_x2(t+h, x1 + k3x1, x2 + k3x2, x3 + k3x3, x4 + k3x4); k4x3 = h * funct_x3(t+h, x1 + k3x1, x2 + k3x2, x3 + k3x3, x4 + k3x4); k4x4 = h * funct_x4(t+h, x1 + k3x1, x2 + k3x2, x3 + k3x3, x4 + k3x4); X1[j+1]= x1+(k1x1 + 2.0 * k2x1 + 2.0 * k3x1 + k4x1) / 6.0; X2[j+1]= x2+(k1x2 + 2.0 * k2x2 + 2.0 * k3x2 + k4x2) / 6.0; X3[j+1]= x3+(k1x3 + 2.0 * k2x3 + 2.0 * k3x3 + k4x3) / 6.0; X4[j+1]= x4+(k1x4 + 2.0 * k2x4 + 2.0 * k3x4 + k4x4) / 6.0; T[j+1] = t + h; j++; } } //---------------------------------------------------------- main () { // matrix for shooting parameters S=(p,q)^T gsl_matrix * S = gsl_matrix_alloc (2, Kmax); // matrix inizializazion for (i = 0; i < 2; i++) for (j = 0; j < Kmax; j++) gsl_matrix_set (S, i, j, 0.0); gsl_vector * s0 = gsl_vector_alloc (2); // initial vector with p and q gsl_vector_set (s0, 0, p); gsl_vector_set (s0, 1, q); gsl_matrix_set_col (S, 0, s0); // matrix for F, F1, F2 gsl_vector * F = gsl_vector_alloc (2); gsl_vector * F1 = gsl_vector_alloc (2); gsl_vector * F2 = gsl_vector_alloc (2); gsl_vector * DF1 = gsl_vector_alloc (2); gsl_vector * DF2 = gsl_vector_alloc (2); gsl_vector * xinv = gsl_vector_alloc (2); gsl_vector * st = gsl_vector_alloc (2); gsl_vector * sM = gsl_vector_alloc (2); gsl_matrix * DF = gsl_matrix_alloc (2, 2); gsl_permutation * p = gsl_permutation_alloc (2); for(k=0; k<=10; k++){ rk4(h, max_steps, 0.0, 1.0, 0.0, gsl_matrix_get (S, 0, k), gsl_matrix_get (S, 1, k)); nn=j; gsl_vector_set (F, 0, X3[nn]-b3); gsl_vector_set (F, 1, X4[nn]-b4); // rk4(h, max_steps, 0.0, 1.0, 0.0, gsl_matrix_get (S, 0, k)+ds, gsl_matrix_get (S, 1, k)); nn=j; gsl_vector_set (F1, 0, X3[nn]-b3); gsl_vector_set (F1, 1, X4[nn]-b4); // rk4(h, max_steps, 0.0, 1.0, 0.0, gsl_matrix_get (S, 0, k), gsl_matrix_get (S, 1, k)+ds); nn=j; gsl_vector_set (F2, 0, X3[nn]-b3); gsl_vector_set (F2, 1, X4[nn]-b4); // Calculate the Jacobi matrix for (i=0; i<=1; i++){ gsl_vector_set (DF1, i, (gsl_vector_get (F1, i)-gsl_vector_get (F, i))/ds); gsl_vector_set (DF2, i, (gsl_vector_get (F2, i)-gsl_vector_get (F, i))/ds); } gsl_matrix_set_col (DF, 0, DF1); gsl_matrix_set_col (DF, 1, DF2); // Calculate the inverse matrix int s; gsl_linalg_LU_decomp (DF, p, &s); gsl_linalg_LU_solve (DF, p, F, xinv); gsl_matrix_get_col (st, S, k); // Calculate the new Newton step gsl_vector_sub (st, xinv); gsl_matrix_set_col (S, k+1, st); } gsl_matrix_get_col (sM, S, k); for(i=0; i<=1; i++){ printf ("The optimal shooting parameters are: S(%d) = %g\n", i, gsl_vector_get(sM, i)); } // Solve equation with the optimal S rk4(h, max_steps, 0.0, 1.0, 0.0, gsl_vector_get (sM, 0), gsl_vector_get (sM, 1)); nn=j; // Plot the resulting solution winsiz (800, 800); page (2600, 2600); sclmod ("full"); scrmod ("revers"); metafl ("xwin"); disini (); pagera (); hwfont (); setrgb(0,0,0); axspos(450,1800); axslen(1800,1200); name("t","x"); name("x(t)","y"); graf(0.0,max_steps,0.0,0.2,0.8,1.1,0.8,0.1); color("blue"); curve(T,X1,nn); endgrf(); // gsl_matrix_free (S); gsl_vector_free (F); gsl_vector_free (F1); gsl_vector_free (F2); gsl_vector_free (s0); gsl_vector_free (DF1);gsl_vector_free (DF2);gsl_matrix_free (DF); disfin (); //---------------------------------------------------------- }