Message Boards Message Boards

0
|
8014 Views
|
4 Replies
|
3 Total Likes
View groups...
Share
Share this post:

Newton-Raphson Root Finding, Difficulty in coding

Posted 12 years ago
Hello, I seem to be having trouble writing for a problem I'm having.  I need to find Specifi Volumes of Carbon Dioxide from the Van der Waals equation, using the Newton method, over a process of constant Temperature and variant Pressure.  I am using Mathematica; I'm not very strong at coding in general and am fairly new to Mathematica.  

I have figured out how to find a single value for the Specific Volume given a set Pressure.  However I am having trouble coding over the pressure change.

I want The pressure to change from 3.9atm to 59atm.  Another issue I seem to be having is it spits our odd answers for Specific Volume if the Temperature is above 329 or the Upper limit of Pressure is above 55.  I need T=333K and the upper Pressure=59.2.  

Here is my code and output, the example pressure is at 10atm:
 a = 1.36;
 b = .003183;
 R = .0820578;
 T = 329;
 P = 10;
 atmA = 3.9;
 atmB = 59.21;
 
 
F[vi_] = vi - {{P + {a/vi^2}} {vi - b} - {R*T}}/{P - {a/vi^2} + {2 a*
       b/vi^3}}
FindRoot[F[vi] == 0, {vi, 5}, WorkingPrecision -> 5]






{{{-((-26.997 + (10 + 1.36/vi^2) (-0.003183 + vi))/(
     10 + 0.00865776/vi^3 - 1.36/vi^2)) + vi}}}

{vi -> 2.6512}

Ignoring these previous issues, my attempt to solve the problem was to use a double integral.  However that has gotten nowehere.
 a = 1.36;
 b = .003183;
 R = .0820578;
 T = 329;
 atmA = 3.9;
 atmB = 59.21;
 
 Integrate[
  vi - {{P + {a/vi^2}} {vi - b} - {R*T}}/{P - {a/vi^2} + {2 a*
       b/vi^3}}, {p, atmA, atmB}, {vi, 0, p}]






{{{4652.6 + 63.0819 I}}}

Thanks for your help.  I know I'm very lost and attempting to teach me may be a hurdle
POSTED BY: Cory Leahy
4 Replies
Posted 12 years ago
I suggest trying the following:
 a = 1.36;
 b = .003183;
 R = .0820578;
 T = 333;
 
 Chop[Array[FindRoot[vi - ((((2.9 + #) + (a/vi^2)) (vi - b) -
   (R T))/((2.9 + #) - (a/vi^2) + (2 a b/vi^3))) == 0,
   {vi, 4 + I}] &, 56]]
 
(* {{vi -> 0}, {vi -> 0.0944562}, {vi -> 0}, {vi -> 0.0944329},
  {vi -> 0}, {vi -> 0.00502882}, {vi -> 0}, {vi -> 0}, {vi -> 0},
  {vi -> 0}, {vi -> 0.0943517}, {vi -> 0.00502902}, ... *)
Please pay attention to the type of parentheses - or braces in your input. Mathematica is very pedantic about it. Although there usually are numerous ways to accomplish a certain goal in Mathematica, input syntax itself doesn't have many synonymous forms.

Second issue that might be relevant here is the fact that FindRoot fails to find real-valued roots for your equation. This is probably because your input parameters are imprecise. Instead of getting bogus results and lots of warnings, I started with complex-valued guess (4 + I), which finds complex-valued roots. Chop removes very small terms from the results - all imaginary parts for the roots are eliminated in this way. I don't know if these results make sense to you; basically root would seem to be more or less precisely zero for every input.

I suggest finding a good tutorial for learning use of Mathematica. It is a bit too large system to learn in a random fashion from the start. Getting a good idea of its' fundamentals takes unfortunately at least couple weeks. I've been personally using it occasional manner - on and off - for almost two decades, and only after couple first years I actually grasped most of the fundamental ideas. Still at this point, I find new stuff at least once a week. It's a powerful system, but like all power tools, tends to demand rigorous approach to learning...
POSTED BY: Jari Kirma
Posted 12 years ago
So I've managed to output a series of equations with Pressures ranging from min to max.  Currently, I'm attempting to take those equations and output a second array which solves for the root, the Specific Volume, of each.  My code:
 a = 1.36;
 b = .003183;
 R = .0820578;
 T = 333;
 
 Array[vi - {{{(2.9 + #) + {a/vi^2}} {vi - b} - {R*T}}/{(2.9 + #) - {a/
          vi^2} + {2 a*b/vi^3}}} &, 56]
 
 Array[FindRoot[% == 0, {vi, 4}], 56]
 Output:
{{{{{-((-27.3252 + (3.9 + 1.36/vi^2) (-0.003183 + vi))/(
       3.9 + 0.00865776/vi^3 - 1.36/vi^2)) +
      vi}}}}, {{{{-((-27.3252 + (4.9 + 1.36/vi^2) (-0.003183 + vi))/(
       4.9 + 0.00865776/vi^3 - 1.36/vi^2)) +
      vi}}}}, {{{{-((-27.3252 + (5.9 + 1.36/vi^2) (-0.003183 + vi))/(
       5.9 + 0.00865776/vi^3 - 1.36/vi^2)) +
      vi}}}}, {{{{-((-27.3252 + (6.9 + 1.36/vi^2) (-0.003183 + vi))/(
       6.9 + 0.00865776/vi^3 - 1.36/vi^2)) + vi}}}},
and the output goes on until the final pressure.

I'm getting an error, "FindRoot::nveq: "The number of equations does not match the number of variables in FindRoot[%==0,{vi,4}" for every equation.  I clearly have the syntax wrong but I don't know how to take the first array and make a second solving for the root of each.
POSTED BY: Cory Leahy
Posted 12 years ago
Jari, thanks for helping.  Upon further searching I believe an Array would be able to satisfy what I need.  I'm thinking I could calculate Specific Volume values over an array of  Van der Waal equations for the given Pressure change.  Another thing I need to do is plot the curve for Pressure vs Specific Volumes.  So this method would work well for getting data points.

However, I've been trying and don't understand the Array format.  For version 8, can someone point me in the right direction for this?
POSTED BY: Cory Leahy
Posted 12 years ago
First of all, try replacing curly braces in your equations (but not in FindRoot or two last arguments of Integrate!) with regular paretheses ( "(" and ")" ). Also, replace = in your definition of F with := .

Curly braces construct lists in Mathematica. Also, when you attempt to define a function, you almost certainly want to use := which is delayed assignment - = is not delayed. Some of the strange output you see is caused by this.

I cannot really comment your actual task in itself, but I suspect something domain-specific is wrong either with your equations, or input parameters. If we do the integration symbolically:
Integrate[vi - ((P + (a/vi^2)) (vi - b) - (R T))/(P - (a/vi^2) + (2 a b/vi^3)),
  {p, a1, a2}, {vi, 0, p},
  Assumptions -> 0 < a1 < a2]

(* ConditionalExpression[(-0.0093055 - 1.35144 a1) a1 + a2 (0.0093055 + 1.35144 a2) + (9.66322*10^-8 - 0.0000151749 a1) Log[0.0063679 - 1. a1] + (0.132705 - 0.363022 a1) Log[0.365557 - 1. a1] + (0.236186 + 0.635037 a1) Log[0.371924 + 1. a1] + (-9.66322*10^-8 + 0.0000151749 a2) Log[0.0063679 - 1. a2] + (-0.132705 + 0.363022 a2) Log[0.365557 - 1. a2] + (-0.236186 - 0.635037 a2) Log[0.371924 + 1. a2], a2 < 0.0063679] *)
The logarithms easily indicate that you are going to get complex-valued result if you input anything higher than 0.0063679 to your pressure range.
POSTED BY: Jari Kirma
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract