Runtime Problem with Simple Diff Eq Solver

[Chris Hammond]

Hi,

I'm using Dev C++ 4.9.9.2 on Windows XP. I'm trying to implement a simple finite difference scheme to solve an ordinary differential equation. The code compiles without issue, but it does not seem to be doing any computation. I create an array to hold the values of my function that I'm trying to find numerically. I input the number of time steps I want to take, and one divided by that is the step-size. With initial conditions u(0)=1 and u'(0)=-1/2, the value of the function at 1 should not be 1, but that's what I'm getting. I don't think anything is happening inside the program. Any help would be appreciated.

Thanks

include <cstdlib>

include <iostream>

include <cmath>

using namespace std;

int main(int argc, char *argv[])
{
float u0=1;
float Du0= -0.5;
int N;
float delta;

cout << "Enter the initial value of u" << endl;
cin >> u0;
cout << "Enter the initial value of u'" << endl;
cin >> Du0;
cout << "Enter the number of steps" << endl;
cin >> N;

delta= 1/N;
float uValues[N+1];
uValues[0]=u0;
uValues[1]=Du0*delta+u0;
cout << uValues[1] << endl;
cout << delta << endl;

for(int count=1; count< N-1;++count){

        uValues[count+1]=(uValues[count]*(2+(1+(count*delta)*(count*delta))
        *delta+(1+count*delta)*delta*delta)+uValues[count-1])/(1+(1+(count*delta)*(count*delta))*delta);

        }
cout << uValues[N];

system("PAUSE");
return EXIT_SUCCESS;

}

[cpns]

It is a matter of type consistency. 1/N = 0 for all N < 1 because both 1 and N have type int. There is an implicit conversion to float only in the assignment after the integer division. So delta becomes zero in all cases.

I suggest that all literal constants that are intended to be floats are declared as such by appending an 'f' and where they are whole numbers (as is the case in all occurances here) '.0f'. I have marked the lines I believe need changing below.

include <cstdlib>

include <iostream>

include <cmath>

using namespace std;

int main(int argc, char *argv[])
{
float u0=1.0f; ////////////////////
float Du0= -0.5f; /////////////////
int N;
float delta;

cout << "Enter the initial value of u" << endl; 
cin >> u0; 
cout << "Enter the initial value of u'" << endl; 
cin >> Du0; 
cout << "Enter the number of steps" << endl; 
cin >> N;

delta= 1.0f / N; //////////////////////// 
float uValues[N+1]; 
uValues[0]=u0; 
uValues[1]=Du0*delta+u0; 
cout << uValues[1] << endl; 
cout << delta << endl;

for(int count=1; count< N-1;++count){

    //////////////////
    uValues[count+1] = (uValues[count]*(2.0f+(1.0f+(count*delta)*(count*delta)) 
        *delta+(1+count*delta)*delta*delta)+uValues[count-1])/(1.0f+(1.0f+(count*delta)*(count*delta))*delta); 
    //////////////////

} 
cout << uValues[N];

system("PAUSE"); 
return EXIT_SUCCESS;

}

I have not tested the code because the line "float uValues[N+1];" is not valid C++ and will not compilein MSVC++ 2008. It is valid in C99 which is not supported by teh MS compiler, and is supported in GCC as an extension to C++. I would avoid such non-portable code is possible. Consider:

float* uValues = new float[N+1];

instead, being sure to later delete it with:

delete [] uValues ;

when done.

I would also consider revising your formatting of the expression in the loop. The appearance of the multiply operator at the beginning of the second line make it look like the unary de-reference operator. As a rule when splitting lines, it is best to leave the binary operator on the preceding line to make it obvious that the line is incomplete and to avoid misinterpretation of the operator (by humans that is - the compiler does not care a jot). Personally I would probably not have produced such a complex expression, and split it into sub expressions assigned to intermediate variables with meaningful names. That approach makes debugging and maintenamce simpler and makes order of evaluation more obvious. Just advice and opinion, you are free to ignore it if you wish (until you ask us to help debug it that is! ;) ).

Clifford

[cpns]

... I made the change to iValues in order to run the test, and got the following output (I have no idea whether the initial values make any sense, or if the result is as you expect)

---

Enter the initial value of u
1
Enter the initial value of u'
2
Enter the number of steps
10
1.2
0.1
-4.31602e+008Press any key to continue . . .

---

Your prompts need some work, and a newline after the result output perhaps.

Clifford

[Chris Hammond]

That was exactly the problem. Thanks. Well, a type-o meant that the answers I got after that were junk, but the real problem was the issue you pointed out.