求助高手Fortran解非线性方程系统

吃草的鸟

五维的方程，xa(5)是要返回的结果，y(5)和其他参数已知，初值赋（/0,0,0,0,0/）
F(1)=xa(1)+y(3)/y(1)+1./cf*(xa(3)-xa(1)**2-1./3.*xa(2)**2-2./3.*xa(2)*y(4)/y(1))
  F(2)=xa(2)+2*y(4)/y(1)+1./cf*(xa(4)-2*xa(1)*xa(2))
  F(3)=xa(3)-1./(1+sigm)*((0.5+0.5*noug/nouf+sigm)*(xa(1)**2+1./3.*xa(2)**2)-noug/nouf*cg*y(3)/y(1)+2./3.*(1+sigm)*xa(2)*y(4)/y(1)-z)
  F(4)=xa(4)-1./(1+sigm)*((1+noug/nouf+2*sigm)*xa(1)*xa(2)-(1-noug/nouf)*g*y(1)-noug/nouf*cg*(2*y(4)/y(1)+y(5)/y(1))+3./8.*Cd*xa(2)**2)
  F(5)=xa(5)-1./(1+sigm)*((2+cg/cf+sigm)*xa(1)*xa(2)+cg/cf*(1-noug/nouf)*g*y(1)-cg*(2*y(4)/y(1)+y(5)/y(1))-cg/cf*3./8.*Cd*xa(2)**2)

我用 call neqnf(fcn,errrel,n,itmax,xguess,x,fnorm)做的，得到的答案和用matlab的fsolve做的不一样啊（matlab的已知是对的）
求：Fortran的正确解方程函数，或者求解方法
大家请赐教！！

吃草的鸟

大家讨论讨论啊，有什么想法说说吧~~

给我提个思路也好，先行谢过！！！

tonnyw

1. Make sure you are using double precision.
2. Use small tolerance.
3. If possible, simplify the equation to be a linear one to see if you get the right  results.

吃草的鸟

3# tonnyw can u speak chinese,please!

吃草的鸟

1. Make sure you are using double precision.
2. Use small tolerance.
3. If possible, simplify the equation to be a linear one to see if you get the right  results.
tonnyw 发表于 2010-8-18 23:05

首先谢谢你的回复

1.double precision ---->赋的值都是real
2.small tolerance   ---->errrel=0.0001
3.非线性怎么转化为线性的啊？

pasuka

首先谢谢你的回复

1.double precision ---->赋的值都是real
2.small tolerance   ---->errrel=0.0001
3.非线性怎么转化为线性的啊？
吃草的鸟 发表于 2010-8-19 11:24

real是单精度
real*8才是双精度！！！

qinxl

是不是Y(1)绝对值太小？如果y(1)等于0，你的方法就会失效。
尽量避免除数中出现接近0的数。

吃草的鸟

是不是Y(1)绝对值太小？如果y(1)等于0，你的方法就会失效。
尽量避免除数中出现接近0的数。
qinxl 发表于 2010-8-19 13:53

xguess（5）赋的初值都是0
初始y=[0.3673,0.0,-68.8183,0,-0.507]

开始matlab和Fortran的结果都是差不多的，后来就不一样了，
解出来后带入一个【龙格库塔4阶】得到结果画图
你帮忙看下这个图，看能发现什么不！

吃草的鸟

real是单精度
real*8才是双精度！！！
pasuka 发表于 2010-8-19 12:57

哦，我看了下书，是单精度的
可是，我已经是保留到小数点后六位了，双精度真的有必要吗？而且errrel=0.0001

qinxl

8# 吃草的鸟

应该跟X(5)的初值关系不大。
大概检查了一下，发现除数中只有y(1),而没有其他4个分量。也就是说，你的算法对于Y是不对称的，好的算法应该不依赖于某一个分量的值。如果你当时凑巧是使用y(2)或者y(4)做除数，那么根据你给的Y值，那么就直接发散了。

建议重组公式，或者重新考虑算法。

吃草的鸟

10# qinxl
谢谢一直以来的指导
还有谢谢大家
我发现程序的错误了，isml的是可以算得，而且我试了试《徐士良》的一个算法，结果是一样的
我想问一下：我定义了有pi=3.14，还有P_I=patm+nouf*9.8*100，为什么循环第一次p_I是对的，第二次以后就是00呢?
我把P_I改为press就好了

qinxl

11# 吃草的鸟
不清楚你的问题啊！

yp51920

很久都没用Fortran了，我觉得你可以从正解去找错误点，很显然正解xa(5)是一个周期循环值，这点很关键；还有把原始的方程给出来给大家看看可能更直观。

吃草的鸟

13# yp51920
方程解出来了
我意思是声明变量是不是不可以使用"_"？？？
比如我定义3.14为pi，压力定义为p_i

yp51920

77中声明变量是不可以使用"_"的
95可以