Community Tip - Learn all about the Community Ranking System, a fun gamification element of the PTC Community. X
Fixed-Point Iteration for nonlinear systems of equations - my Mathcad alghoritm doesn't work
Hello,
I try to make my Mathcad alghoritm to work for a nonlinear system by using Fixed-Point Iteration . Unfortunately still I cannot get the results. Something is missing in my script. I attached my Mathcad file.
The results whould look finally like this (I obtained it in Excel in order to compare the results if possible):
It xold yold x y errx erry errf10 errf20
-----------------------------------------------
1 2.0000 4.0000 1.4142 6.0828 0.5858 2.0828 0.6023 23.7569
2 1.4142 6.0828 1.6132 7.3348 0.1990 1.2520 3.2299 9.4517
3 1.6132 7.3348 2.4150 6.9238 0.8018 0.4109 4.8888 38.8502
4 2.4150 6.9238 3.2743 5.6438 0.8593 1.2800 1.7585 28.6511
5 3.2743 5.6438 3.5326 4.8065 0.2583 0.8373 1.5000 6.8054
6 3.5326 4.8065 3.3135 4.6017 0.2191 0.2048 1.7316 4.4348
7 3.3135 4.6017 3.0410 4.7449 0.2725 0.1432 0.8185 6.2781
8 3.0410 4.7449 2.9034 4.9577 0.1377 0.2128 0.0356 3.5969
9 2.9034 4.9577 2.8972 5.0811 0.0061 0.1234 0.3272 0.2816
10 2.8972 5.0811 2.9532 5.0907 0.0559 0.0096 0.3124 1.4396
11 2.9532 5.0907 3.0056 5.0425 0.0524 0.0481 0.1223 1.3810
12 3.0056 5.0425 3.0258 4.9968 0.0203 0.0458 0.0363 0.5519
13 3.0258 4.9968 3.0198 4.9785 0.0060 0.0183 0.0852 0.1304
14 3.0198 4.9785 3.0057 4.9828 0.0141 0.0043 0.0574 0.3554
15 3.0057 4.9828 2.9962 4.9947 0.0096 0.0119 0.0121 0.2503
16 2.9962 4.9947 2.9941 5.0030 0.0020 0.0084 0.0149 0.0590
17 2.9941 5.0030 2.9966 5.0050 0.0025 0.0020 0.0184 0.0604
18 2.9966 5.0050 2.9997 5.0030 0.0031 0.0020 0.0093 0.0786
19 2.9997 5.0030 3.0012 5.0004 0.0015 0.0026 0.0001 0.0413
20 3.0012 5.0004 3.0012 4.9990 0.0000 0.0014 0.0042 0.0008
21 3.0012 4.9990 3.0005 4.9990 0.0007 0.0000 0.0036 0.0176
Now, the code is rigorous also for programers. I initialized k1,k2, k3,k4 -> 0.
Any of the ways is chosen has its own inconvenients, by mathematical, programming or of mode of displaying point of view. The programming environment is always a limited one.
And for Midpoint Method, which one do you prefer Werner ?
I tend to think that the one from the left is simpler. Still it requires the trick to put y after the Tab reading data. Cool I learnt several programming tricks from you last days.
Yes, I definitely prefer the left variant. Not only because its a bit cleaner and simpler, but also because I like that the meaning of the function arguments remains the same throughout the whole program.
In your right viriant x0 and y0 initially have the meaning of start value( initial condition. But later x0 and y0 the meaning changes and they are used for something like next_x, next_y.
I still don't understand why you insists on using a and b in your program but do not want to name the function arguments accordingly.
Also the assignment x<--x0 could be deleted - its not necessary as x is calculated as x0+h*i anyway. This saves two more program lines. Don't get me wrong - I don't think that a program is better just because its shorter. A program with more lines which is better readable and understandable often is preferrable. But I guess this does not apply to your right variant Midpoint1.
My personal preference is to provide the number of steps as last argument and not the step width. This ensures that we really reach the end value. Try a value of h=0.08 in one of your variants and see what I mean.
Given that, the program could be shortened once again 😉
The reason I do not calculate h and use it in the loop is to avoid summation of round off errors. The calculation is so simple and fast that I feel its reasonable to let the program do the calculation many times instead of just once.
And for Heun method can be added somehow inside FOR i=0..n instruction a specification with IF i=0, x ->x0, y-> y0, else - > the main algorithm ? This in order to include the initial conditions in the table (without vectorization).
I found something, but the problem remains if I can simplify more the code. I want to use more x, y, yp instead of terms like x0,y0, etc.
Something like this
MC15 sheet attached
This is the simplest variant i am working on but still don't work. It eliminates several x0, y0. The scheme runs OK in C++, Matlab but I still cannot make it work in Mathcad. I am sure it is good, but Mathcad is still opposing.
I get closing on the solution... But the solution is still a bit wrong.
I found ! There are stages where are necessary 2 types of x from the same iteration. An old and an actual x. Now it is working with the good solution.
@tubar wrote:
I found ! There are stages where are necessary 2 types of x from the same iteration. An old and an actual x. Now it is working with the good solution.
Hmm, so all you were looking after was that you wanted f(x,y) to show the value from x and y of the same row, but yp and f(x,yp) should be derived from the values from the row above it? Seems not logical to me but if thats what you want, you obviously finally found a solution.
If you really like that output, here is simpler way to achieve it:
IMHO it would make more sense to put the "yp" and "f(x,yp)" columns in front before the "x" and "y" columns.
@tubar wrote:
This is the simplest variant i am working on but still don't work. It eliminates several x0, y0. The scheme runs OK in C++, Matlab but I still cannot make it work in Mathcad. I am sure it is good, but Mathcad is still opposing.
The red error tells it - you are using yp but that variable never is assigned a value, so its unknown. This sure would yield a similar error in C++ or Matlab, too.
For trapezoid integration, do you know why the solution from the right doesn't work ? It implements my scheme philosophy from last days.
I found now this way. I was forgetting to respect an order of the equations.
I try now to implement to all my files the philosophy from last days.
With "philosophy" you mean the if..otherwise... statements, I guess.
I see no reason for the first for-loop in the left variant and a see no reason for making x and y vectors in both variants.
Finally I found the solution. I put it even to find out the area under the graph.
I split FOR in two stages based on IF.
It is interesting how you escaped from vectors, because I tried a lot.
Looks quite elaborate.
According avoiding the vectors - The only value you want to keep because your algorithm is using it multiple times (could also be avoided) is y0. It didn't seem reasonable to me to create vectors for x and y just to have access to y0. So I made y0 a separate variable.
I finally succeed for Simpson. An d I will try also for the other Simpson 3-8. It evaluates also the intermediary cumulative values.
It is a way to put the IF condition not on the same line with the consequence equation when that equation is only one? I want to save some horizontal space.
IF (condition)
| Consequence equation
instead of :
Consequence equation IF (condition)
@tubar wrote:
It is a way to put the IF condition not on the same line with the consequence equation when that equation is only one? I want to save some horizontal space.
IF (condition)
| Consequence equation
instead of :
Consequence equation IF (condition)
There is no way to override Mathcads autofomatting of the if-command, but you can make a two-line consequence by adding a dummy command, just a number or even better just add a line with a string (either a null-string or use it for a comment) so you have a two-line "consequence.
The one line if-command saves vertical space, but I, too, find it often more important to save horizontal space and I think that the multiple line variant provides better readability, too.
Thanks. Werner, I want to invite you in the future to collaborate to some of my research articles (domain engineering thermodynamics) if you transmit to me your email address at tbaracu@yahoo.co.uk .
There may be some complex equations solvable through Mathcad where you can be useful.
My knowledge in the field you are working in is close to zero, but if you get stuck with the realization of a problem in Mathcad, feel free to post your question in the forum here. Of course it may help if you don't presuppose knowledge in thermodynamics and try to formulate the question with the focus on the pure math 😉
for this Gauss Quadrature integration, do you identify the origin of the error, Werner ? It is somehow in a block that I named Correctness, but I dont know why it cannot apply Newton-Raphson nonlin. eq. solving.
I attached the Mathcad file and also the C++/txt file for good orientation.
The final solution shoud appear like this (obtained in C++):
Enter the value of n for Pn(x) :
4
Enter the lower limit a of integration :
0
Enter the upper limit b of integration :
1
The Legendre's Polynomial is : 0.375 + (-3.75) X^2 + (4.375) X^4
Roots x(i) Weights w(i)
0.861136311594053 0.347854845137454
0.339981043584856 0.652145154862546
-0.339981043584856 0.652145154862546
-0.861136311594053 0.347854845137454
The Value of Integration is I= 0.132120558485338
Haven't looked for the cause but the error occurs in your function Pn which is called by CheckCorrectness with a value of x in the range of 10^91.
You may turn debug mode on and throw in some trace commands in your functions to track down why x is growing up to infinite.
It is necessary inside the sub-block CheckCorrectness a kind of vectorialization (through i) inside the While Loop that is inside the FOR loop. C++ tolerated simple variables, but Mathcad may require this vectorialization.
I will try to repair initially the sub-block CheckCorectness be taking it separately, isolating it, considering Pn(n,x) and its derivate Pndx(n,x) with their formula (taken with an arbitrary form, as example) exp(-x)-x and -exp(-x)-1 respectively.
This is the isolated sub-block Checkcorrectness (attached). I try to solve it separately as a Newton-Raphson method.
I solved that sub-block Correctness. It needed an intermediary variable inside FOR loop to read the values of x. I will reconsider now the main file.
I think is a problem with Mathcad internal convergence, 10^-16 is too much for it. I try to set it other value. How you obtained DEBUG info ? I presed Toggle Info but it shew nothing.
When you press "Toggle Debug" you should get the Trace Window at the bottom of the screen.
You have to throw in the "trace" command in your program (look it up in the help) and then let the sheet recalculate (Crtl-F9)
And yes, 10^-16 is too much precision demanded for a 32 bit program with variables in IEEE format. This also applies to a C++ program (at least when using a 32bit compiler).
Yeeeeeeessss ! I finally obtained the results in the main file !!!!!! It is able to calculate itself the POLINOMIALS, the ROOTS and the WEIGHTING FACTORS. You can choose any order of Gauss Quadrature.
I used also the SYMBOLIC BLOCK to verify the formula of Pn, Pndx .
So It is no way to make Mathcad allow max errors of 10^16 . When i set 10^-15 it works. When I set it 10^-16 it is processing in idle mode.
In C++ not problem ! I run the file with DevC++ and also Codeblocks. It just did the job.