0

I´ve been using this code in C for integration with different equations but today i modified it to integrate the one shown and the .dat file is giving me all the columns full of "-nan". Is it a wrong-coding issue or its just that this equations isnt meant to be solved by this procedure?

Here are both the integration routine and the main code. Thanks

FIRST THE CODE (i chopped it explaining what does each part do, i apologize if its unpleasant at sight)

#include <stdio.h> #include <math.h> #define a -1

PARAMETERS

struct Par{ double mu1, mu2, w1, w2, eps; } aa;

EQUATIONS

void ecuaciones(int n, double v[], double dv[], double t){ double x1,y1,x2,y2; x1=v[0]; x2=v[1]; y1=v[2]; y2=v[3]; y1=aa.mu1*x1-aa.w1*y1-(x1*x1+y1*y1)*x1+aa.eps*(x1-x2) ; y2=aa.mu2*x2-aa.w2*y2-(x2*x2+y2*y2)*x2+aa.eps*(x2-x1) ; dv[0]=y1; dv[1]=y2; dv[2]=aa.w1*x1+aa.mu1*y1-y1*(x1*x1+y1*y1)+aa.eps*(y1-y2) ; dv[3]=aa.w2*x2+aa.mu2*y2-y2*(x2*x2+y2*y2)+aa.eps*(y2-y1) ; return; }

MAIN

int main(){ int i,j; FILE *ptr; double v[4],t,dt,t_pre,t_max;

EXIT

 ptr=fopen("NLIAC.dat","w"); dt=0.01; t_max=20;

INITIAL CONDITIONS

for (i = 2; i < 6; i++) { v[0]=i; v[1]=6; v[2]=0; v[3]=1;

PARAMETERS DEFINITION aa.w1=1; aa.mu1=1; aa.w2=6; aa.mu2=1; aa.eps=0;

 t=0.;

INTEGRATION COMMAND

while(t<t_max){ rk4(ecuaciones,v,4,t,dt);

EXPORT

fprintf(ptr,"%lg\t%lg\t%lg\t%lg\t%lg\n",t,v[0],v[1],v[2],v[3]); t+=dt; }} fprintf(ptr,"\n"); fclose(ptr); return(0); }

AND HERE IS THE INTEGRATION ROUTINE (its fine)

 void rk4(void deri(int , double [], double [], double ), \ double h[], int n, double t, double dt) { #define naux 26 int i; double k1[naux],k2[naux],k3[naux],k4[naux],h0[naux]; double dt2, dt6; dt2=dt/2.; dt6=dt/6.; for (i = 0 ; i<n; i++) h0[i] = h[i]; deri(n,h0,k1,t); for (i =0 ; i<n ; i++) h0[i]=h[i]+dt2*k1[i]; deri(n,h0,k2,t+dt2); for (i =0 ; i<n ; i++) h0[i]=h[i]+dt2*k2[i]; deri(n,h0,k3,t+dt2); for (i =0 ; i<n ; i++) h0[i]=h[i]+dt*k3[i]; deri(n,h0,k4,t+dt); for (i = 0 ; i<n ; i++) {h0[i]=h[i]+dt*k4[i];}; for (i =0; i<n ; i++) h[i]=h[i]+dt6*(2.*(k2[i]+k3[i])+k1[i]+k4[i]); return; }
Improve this question
4
  • I would check on the re-definition of y1 and y2. You are using y1 in the second formula, however its value has changed. Define dv[0] and dv[1] directly without the intermediate assignment. Commented Jul 13, 2016 at 1:51
  • thanks i will try and comment how it turned out Commented Jul 13, 2016 at 2:02
  • I am pleased to announce that what LutzL suggested is indeed the way to solve the problem. This thread can be closed. Thanks for your help!!! Commented Jul 13, 2016 at 2:15
  • @GatoMazzei couldn't you have posted that 2 minutes earlier? ;-) But serious: if LutzL does not want to post an answer you can answer your question yourself. Commented Jul 13, 2016 at 2:23

1 Answer 1

Reset to default
1

You get a massive overflow in y2 quite quickly which avalanches even faster into NaNs

Never undestimate the usefulness of printfs in numerical code, so let's do it here:

void ecuaciones(int n, double v[], double dv[], double t) { double x1, y1, x2, y2; x1 = v[0]; x2 = v[1]; y1 = v[2]; y2 = v[3]; printf("aa.mu1: %g, aa.w1: %g, aa.eps: %g, aa.mu2: %g, aa.w2: %g \n", aa.mu1, aa.w1, aa.eps, aa.mu2, aa.w2); printf("x1: %g, x2: %g, y1: %g, y2: %g\n", x1, x2, y1, y2); y1 = aa.mu1 * x1 - aa.w1 * y1 - (x1 * x1 + y1 * y1) * x1 + aa.eps * (x1 - x2); y2 = aa.mu2 * x2 - aa.w2 * y2 - (x2 * x2 + y2 * y2) * x2 + aa.eps * (x2 - x1); dv[0] = y1; printf("y1: %g, y2: %g\n", y1, y2); dv[1] = y2; dv[2] = aa.w1 * x1 + aa.mu1 * y1 - y1 * (x1 * x1 + y1 * y1) + aa.eps * (y1 - y2); dv[3] = aa.w2 * x2 + aa.mu2 * y2 - y2 * (x2 * x2 + y2 * y2) + aa.eps * (y2 - y1); printf("%g %g %g %g %g\n\n", dv[0], dv[1], dv[2], dv[3], y1); return; }

shows:

aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 x1: 2, x2: 6, y1: 0, y2: 1 y1: -6, y2: -222 -6 -222 236 1.09489e+07 -6 aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 x1: 1.97, x2: 4.89, y1: 1.18, y2: 54745.3 y1: -9.5984, y2: -1.46559e+10 -9.5984 -1.46559e+10 913.916 3.148e+30 -9.5984 aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 x1: 1.95201, x2: -7.32794e+07, y1: 4.56958, y2: 1.574e+28 y1: -50.8154, y2: 1.81548e+64 -50.8154 1.81548e+64 131360 -5.98381e+192 -50.8154 aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 x1: 1.49185, x2: 1.81548e+62, y1: 1313.6, y2: -5.98381e+190 y1: -2.57558e+06, y2: -inf -2.57558e+06 -inf -nan -nan -2.57558e+06 aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 x1: -4290.84, x2: -inf, y1: -nan, y2: -nan y1: -nan, y2: -nan -nan -nan -nan -nan -nan aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 x1: -nan, x2: -nan, y1: -nan, y2: -nan y1: -nan, y2: -nan -nan -nan -nan -nan -nan
Improve this answer
Sign up to request clarification or add additional context in comments.

6 Comments

Great answer! not only explaining what is wrong but teaching what a great tool printf is!
Now: is it valid to use "v[number]" as a variable in the equation (dv=..), that way avoiding defining a name for each one?
Why do you expect to get an answer after making fun of my? ;-) You can, of course, use the elements of the vector directly, no need for a temporary variable. Using a temporary variable also lead me to the wrong conclusion that all of the overwriting of yn was intentionally and the overflow was the problem instead of "just" a typo. And, yes, I am serious with the printfs in numerical computations, you can even plot them, given enough data points, and find e.g.: a "nervous converging" (for lack of a better word) which is hard to find otherwise.
Actually im not a huge fan of sarcasm. That means i meant exactly what i said! :). The remaining thing in my work is now determine if the -nan obtained for certain initial conditions are part of the physical phenomena being studied or there is something else wrong with the code. Given the fact that this is supposed to be my final exam for tomorrow i´ll probably ignore this and present an average-joe work :(. The two lessons that readers can learn from this post are manily:
1) start your numerical simulations as soon as possible: you never know what problems will occur
|

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.