2019年2月28日星期四

Julia語言習題3-6請選擇適當的微分近似法完成下表：

Julia語言習題3-6請選擇適當的微分近似法完成下表：
x = 0 ,2 ,3 ,5
f(x)= 1 ,7 , 25 ,121
f'(x) = ? ? ? ?

程式
using Printf

x=[ 0 , 2 , 3 , 5]
f= [ [1.0 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[7.0 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[25.0 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[121.0 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ] ]

n=length(x)
println(" Divided Difference Table: ")
println("=============================")
for j=2:n
for i=1:n-j+1
f[i][j]=(f[i+1][j-1]-f[i][j-1])/(x[i+j-1]-x[i])
end
end

print("i\tx(i)\t\tf(i)\t\tf(i,i+1)\tf(i,i+1.i+2), ......................\n")
for i=1:n
s=@sprintf("%d\t%8.5f",i,x[i])
print(s)
for j=1:n-i+1
s=@sprintf("\t%8.5f",f[i][j])
print(s)
end
println()

end

d=[0.0 for i=1:n-1]
for i=1:n-1
d[i]=f[1][i+1]
end

s=@sprintf("\n牛頓一階微分前向差除表")
println(s,d)

xa=[ 0.0 , 2.0 , 3.0 , 5.0]
pn=[0.0 , 0.0 , 0.0 ,0.0]
s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
pn[i]=d[1]
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end

s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
xb=xa[i]
dx=0.0
for j=1:2
#println(xb," ", x[j])
dx=dx+(xb-x[j])
end
#println(dx)
dx=d[2]*dx
pn[i]=pn[i]+dx
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end

s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
xb=xa[i]
dx=0.0
dx=(xb-x[1])*(xb-x[2])+ (xb-x[2])*(xb-x[3])+ (xb-x[3])*(xb-x[1])
dx=d[3]*dx
pn[i]=pn[i]+dx
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end

輸出畫面

Divided Difference Table: 
=============================
i x(i)  f(i)  f(i,i+1) f(i,i+1.i+2),  ......................
1  0.00000  1.00000  3.00000  5.00000  1.00000
2  2.00000  7.00000 18.00000 10.00000
3  3.00000 25.00000 48.00000
4  5.00000 121.00000

牛頓一階微分前向差除表[3.0, 5.0, 1.0]

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.000 , P'n( 0.000 ) = 3.00000
========================================================
xa[2]= 2.000 , P'n( 2.000 ) = 3.00000
========================================================
xa[3]= 3.000 , P'n( 3.000 ) = 3.00000
========================================================
xa[4]= 5.000 , P'n( 5.000 ) = 3.00000

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.000 , P'n( 0.000 ) = -7.00000
========================================================
xa[2]= 2.000 , P'n( 2.000 ) = 13.00000
========================================================
xa[3]= 3.000 , P'n( 3.000 ) = 23.00000
========================================================
xa[4]= 5.000 , P'n( 5.000 ) = 43.00000

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.000 , P'n( 0.000 ) = -1.00000
========================================================
xa[2]= 2.000 , P'n( 2.000 ) = 11.00000
========================================================
xa[3]= 3.000 , P'n( 3.000 ) = 26.00000
========================================================
xa[4]= 5.000 , P'n( 5.000 ) = 74.00000

Julia語言微分的近似值習題3-3 完成表中的 f'(x) 與 f''(x)

Julia語言微分的近似值習題3-3 完成表中的 f'(x) 與 f''(x)
微分的近似值完成表中的 f'(x) 與 f''(x)
x= f(x)= f'(x)= f''(x)=
=============================
0.15 0.1761 ?? ??
0.17 0.2304 ?? ??
0.19 0.2788 ?? ??
0.21 0.3222 ?? ??
=============================
h=0.02

using Printf

x=[ 0.15 , 0.17 , 0.19 , 0.21 ]
f= [ [0.1761 ,0.0 , 0.0 ],
[0.2304 ,0.0 , 0.0 ],
[0.2788 ,0.0 , 0.0 ],
[0.3222 ,0.0 , 0.0 ] ]

h=abs(x[1]-x[2])
i=1
s=@sprintf("f(%0.3f)的微分.......",x[i])
println(s)

fd1=( 2*f[i+3][1] - 9*f[i+2][1] + 18*f[i+1][1] - 11*f[i][1] ) / (6*h)
fd2=( - f[i+3][1] + 4*f[i+2][1] + -5*f[i+1][1] + 2*f[i][1] ) / (h*h)
s=@sprintf("一次微分值= %0.4f , 二次微分值= %0.4f ",fd1,fd2)
f[i][2]=fd1
f[i][3]=fd2
println(s)

println("\n")
i=2
s=@sprintf("f(%0.3f)的微分.......",x[i])
println(s)

fd1=(f[i+1][1]- f[i+-1][1]) / (2*h)
fd2=(f[i+1][1] - 2*f[i][1] + f[i-1][1] ) / (h*h)
s=@sprintf("一次微分值= %0.4f , 二次微分值= %0.4f ",fd1,fd2)
f[i][2]=fd1
f[i][3]=fd2
println(s)

println("\n")
i=3
s=@sprintf("f(%0.3f)的微分.......",x[i])
println(s)
fd1=(f[i+1][1]- f[i+-1][1]) / (2*h)
fd2=(f[i+1][1] - 2*f[i][1] + f[i-1][1] ) / (h*h)
s=@sprintf("一次微分值= %0.4f , 二次微分值= %0.4f ",fd1,fd2)
f[i][2]=fd1
f[i][3]=fd2
println(s)

println("\n")
i=4
s=@sprintf("f(%0.3f)的微分.......",x[i])
println(s)

fd1=( 11*f[i][1] - 18*f[i-1][1] + 9*f[i-2][1] - 2*f[i-3][1] ) / (6*h)
fd2=( 2*f[i][1] - 5*f[i-1][1] + 4*f[i-2][1] - f[i-3][1] ) / (h*h)
s=@sprintf("一次微分值= %0.4f , 二次微分值= %0.4f ",fd1,fd2)
f[i][2]=fd1
f[i][3]=fd2
println(s)

println("\n")
println("x\tf(x)\t\tf'(x)\t\tf''(x)")
for i=1:4
print(x[i],"\t")
s=@sprintf("%0.5f\t\t%0.5f\t\t%0.5f",f[i][1],f[i][2],f[i][3])
println(s)
end

輸出畫面

f(0.150)的微分.......
一次微分值= 2.8775 , 二次微分值= -17.0000 


f(0.170)的微分.......
一次微分值= 2.5675 , 二次微分值= -14.7500 


f(0.190)的微分.......
一次微分值= 2.2950 , 二次微分值= -12.5000 


f(0.210)的微分.......
一次微分值= 2.0600 , 二次微分值= -10.2500

x     f(x)     f'(x)     f''(x)

0.15 0.17610  2.87750  -17.00000
0.17 0.23040  2.56750  -14.75000
0.19 0.27880  2.29500  -12.50000
0.21 0.32220  2.06000  -10.25000

Julia語言例題3-3 不等距的函數f(x)的微分近似值先使用牛頓向前的內插多項式

Julia語言例題3-3 不等距的函數f(x)的微分近似值
先使用牛頓向前的內插多項式

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝====
x f(x) 第一項P'n(x) 前二項P'n(x)　　前三項P'n(x)
0.5 0.4794
0.6 0.5646
0.8 0.7174
1.05 0.8674
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝====
P'(n) 第1項 f0,1

P'(n) 第2項 f0,1+f0,1,2 *( (x-x[1]) + (x-x[2]) )

P'(n) 第3項 f0,1+f0,1,2 *( (x-x[1]) + (x-x[2]) )
+ f0,1,2,3 * ( (x-x[1])*(x-x[2]) ) + ( (x-x[2])*(x-x[3]) ) + ( (x-x[3])*(x-x[1]) )

using Printf

x=[ 0.5 , 0.6 , 0.8 , 1.05]
f= [ [0.4794 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[0.5646 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[0.7174 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[0.8674 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ] ]

n=length(x)
println(" Divided Difference Table: ")
println("=============================")
for j=2:n
for i=1:n-j+1
f[i][j]=(f[i+1][j-1]-f[i][j-1])/(x[i+j-1]-x[i])
end
end

print("i\tx(i)\t\tf(i)\t\tf(i,i+1)\tf(i,i+1.i+2), ......................\n")
for i=1:n
s=@sprintf("%d\t%8.5f",i,x[i])
print(s)
for j=1:n-i+1
s=@sprintf("\t%8.5f",f[i][j])
print(s)
end
println()

end

d=[0.0 for i=1:n-1]
for i=1:n-1
d[i]=f[1][i+1]
end

s=@sprintf("\n牛頓一階微分前向差除表")
println(s,d)

xa=[ 0.5 , 0.6 , 0.8 , 1.05]
pn=[0.0 , 0.0 , 0.0 ,0.0]
s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
pn[i]=d[1]
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end

s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
xb=xa[i]
dx=0.0
for j=1:2
#println(xb," ", x[j])
dx=dx+(xb-x[j])
end
#println(dx)
dx=d[2]*dx
pn[i]=pn[i]+dx
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end

s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
xb=xa[i]
dx=0.0
dx=(xb-x[1])*(xb-x[2])+ (xb-x[2])*(xb-x[3])+ (xb-x[3])*(xb-x[1])
dx=d[3]*dx
pn[i]=pn[i]+dx
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end

println("\n\n")
for i=1:length(xa)
xb=xa[i]
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , 誤差 f'(x)-P'(n) = %0.5f",i,xa[i],abs(cos(xb)-pn[i]) )
println(s)
end

輸出畫面

i f(i)  f(i,i+1) f(i,i+1.i+2),  ......................
1  0.50000  0.47940  0.85200 -0.29333 -0.12929
2  0.60000  0.56460  0.76400 -0.36444
3  0.80000  0.71740  0.60000
4  1.05000  0.86740

牛頓一階微分前向差除表[0.852, -0.293333, -0.129293]

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.500 , P'n( 0.500 ) = 0.85200
========================================================
xa[2]= 0.600 , P'n( 0.600 ) = 0.85200
========================================================
xa[3]= 0.800 , P'n( 0.800 ) = 0.85200
========================================================
xa[4]= 1.050 , P'n( 1.050 ) = 0.85200

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.500 , P'n( 0.500 ) = 0.88133
========================================================
xa[2]= 0.600 , P'n( 0.600 ) = 0.82267
========================================================
xa[3]= 0.800 , P'n( 0.800 ) = 0.70533
========================================================
xa[4]= 1.050 , P'n( 1.050 ) = 0.55867

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.500 , P'n( 0.500 ) = 0.87745
========================================================
xa[2]= 0.600 , P'n( 0.600 ) = 0.82525
========================================================
xa[3]= 0.800 , P'n( 0.800 ) = 0.69758
========================================================
xa[4]= 1.050 , P'n( 1.050 ) = 0.49434

========================================================
xa[1]= 0.500 , 誤差 f'(x)-P'(n) = 0.00013
========================================================
xa[2]= 0.600 , 誤差 f'(x)-P'(n) = 0.00008
========================================================
xa[3]= 0.800 , 誤差 f'(x)-P'(n) = 0.00087
========================================================
xa[4]= 1.050 , 誤差 f'(x)-P'(n) = 0.00323

2019年2月25日星期一

Julia語言利用3點微分方式求f'(x)

Julia語言利用3點微分方式求f'(x)
3點向前
fd1=((-f[i+2]+ 4*f[i+1]-3*f[i])/ (2*h))
3點向後
fd2=((3*f[m] -4*f[m-1]+ f[m-2]) / (2*h))

x f(x) f'(x)
======================
-0.3 -0.20431 ???
-0.1 -0.08993 ???
0.1 0.11007 ???
0.3 0.39569 ???
======================

using Printf

#=====================================================
Trigonometric and hyperbolic functions
All the standard trigonometric and hyperbolic functions are also defined:

sin cos tan cot sec csc
sinh cosh tanh coth sech csch
asin acos atan acot asec acsc
asinh acosh atanh acoth asech acsch
sinc cosc atan2
=====================================================#

x=[-0.3 , -0.1, 0.1 , 0.3 ]
f=[-0.20431 , -0.08993 , 0.11007 , 0.39569]

println("3點向前差")
for i=1: (length(x)-2)
h= abs(x[i+1]-x[i])
fd1=((-f[i+2]+ 4*f[i+1]-3*f[i])/ (2*h))
s=@sprintf("3點向前差 Three-point fordward Diff : fd(%0.5lf) =%0.5lf \n",x[i], fd1)
println(s)
end

println("3點向後差")
for m=length(x):-1:3
h= abs(x[m]-x[m-1])
fd2=((3*f[m] -4*f[m-1]+ f[m-2]) / (2*h))
s=@sprintf("3點向後差 Three-point bachfoward Diff : fd(%0.5lf) =%0.5lf \n",x[m], fd2)
println(s)
end

輸出畫面

3點向前差
3點向前差 Three-point fordward Diff : fd(-0.30000) =0.35785 

3點向前差 Three-point fordward Diff : fd(-0.10000) =0.78595 

3點向後差
3點向後差 Three-point bachfoward Diff : fd(0.30000) =1.64215 

3點向後差 Three-point bachfoward Diff : fd(0.10000) =1.21405

Julia語言例題3-1 利用向前差　向後差　中央差　求f(x)=sin(x)的值

Julia語言例題3-1 利用向前差　向後差　中央差　求f(x)=sin(x)的值
h=0.001 , 0.005 , 0.01 , 0.05 , 0.1 ,0.5

def func( x): # // 欲微分函數
return math.sin(x)

def FordDiff( x, h, fx):
# // 前差微分
return ( ( fx(x+h) - fx(x) ) / h)

def BackDiff(x, h, fx):
# // 後差微分
return (( fx(x) - fx(x-h) ) / h)

def MidDiff(x, h, fx):
# // 中差微分
return (0.5 * ( fx(x+h) - fx(x-h) ) / h)

using Printf

#=====================================================
Trigonometric and hyperbolic functions
All the standard trigonometric and hyperbolic functions are also defined:

sin cos tan cot sec csc
sinh cosh tanh coth sech csch
asin acos atan acot asec acsc
asinh acosh atanh acoth asech acsch
sinc cosc atan2
=====================================================#

function func(x::Float64) #// 欲微分函數
return sin(x)
end

function FordDiff(x::Float64 , h::Float64, fx::typeof(func)) #//前差微分
a=fx(x+h)
b=fx(x)
return ((a-b) / h)
end

function BackDiff(x::Float64 , h::Float64, fx::typeof(func)) #//後差微分
return ((fx(x) - fx(x-h) ) / h )
end

function MidDiff(x::Float64 , h::Float64, fx::typeof(func)) #//中差微分
return (0.5 * ( fx(x+h) - fx(x-h) ) / h)
end

#=====================================================#
hx=[0.001 , 0.005 , 0.01 , 0.05 , 0.1 , 0.5]

for m=1:length(hx)
x=1.0
h=hx[m]
println("h=",h,"\n")
answer = cos(x) # // 答案
cal = FordDiff(x, h, func)
delta = (cal - answer)/answer
s=@sprintf("FordDiff : %0.5lf, delta = %0.5lf %%\n",cal,delta)
println(s)

cal = BackDiff(x, h, func)
delta = (cal - answer)/answer
s=@sprintf("BackDiff : %0.5lf, delta = %0.5lf %%\n",cal,delta)
println(s)

cal = MidDiff(x, h, func)
delta = (cal - answer)/answer
s=@sprintf("MidDiff : %0.5lf, delta = %0.5lf %%\n",cal,delta)
println(s)
println("================================================")
end

輸出畫面

h=0.001

FordDiff : 0.53988, delta = -0.00078 %

BackDiff : 0.54072, delta = 0.00078 %

MidDiff : 0.54030, delta = -0.00000 %

================================================
h=0.005

FordDiff : 0.53820, delta = -0.00390 %

BackDiff : 0.54240, delta = 0.00389 %

MidDiff : 0.54030, delta = -0.00000 %

================================================
h=0.01

FordDiff : 0.53609, delta = -0.00780 %

BackDiff : 0.54450, delta = 0.00777 %

MidDiff : 0.54029, delta = -0.00002 %

================================================
h=0.05

FordDiff : 0.51904, delta = -0.03934 %

BackDiff : 0.56111, delta = 0.03851 %

MidDiff : 0.54008, delta = -0.00042 %

================================================
h=0.1

FordDiff : 0.49736, delta = -0.07947 %

BackDiff : 0.58144, delta = 0.07614 %

MidDiff : 0.53940, delta = -0.00167 %

================================================
h=0.5

FordDiff : 0.31205, delta = -0.42246 %

BackDiff : 0.72409, delta = 0.34016 %

MidDiff : 0.51807, delta = -0.04115 %

================================================

Julia語言例題 EX2-8 定點回路法求非線性方程式 f(x)=exp(x) - 3x^2 = 0

Julia語言例題 EX2-8 定點回路法求非線性方程式 f(x)=exp(x) - 3x^2 = 0

'''
定點回路法
/* ex2-8.jl is used for solving nonlinear equation
* based on Fixed-Point Algorithm g(x)=x with initial
* approximation P0.
*/

定點回路法求非線性方程式 f(x)=math.exp(x) - 3x^2 = 0
可以改寫成
y= x
y= (esp(x)/3) ^ 0.5

與
y= x
y= - (esp(x)/3) ^ 0.5

using Printf

#=======================================================
* ex2-8.jl is used for solving nonlinear equation
* based on Fixed-Point Algorithm g(x)=x with initial
* approximation P0.
=========================================================#
MAX=50
TOL=0.0001

function g1(x0::Float64)
return ( (exp(x)/3)^0.5)
end

function g2(x0::Float64)
return (-(exp(x)/3)^0.5)
end

i=1
x0=0.0
x=0.0
while(i<=MAX)
x=g1(x0)
s=@sprintf("%-2d %10.7lf",i-1,x0)
println(s)
if(abs(x-x0) < TOL)
s=@sprintf("The Root=%10.7lf x-x0=%10.7lf",x,abs(x-x0))
println(s)
break
end
i+=1
x0=x
end
s=@sprintf("Fixed-point failed after %d iteration.\n",i)
println(s)

i=1
x0=0.0
x=0.0
while(i<=MAX)
x=g2(x0)
s=@sprintf("%-2d %10.7lf",i-1,x0)
println(s)
if(abs(x-x0) < TOL)
s=@sprintf("The Root=%10.7lf x-x0=%10.7lf",x,abs(x-x0))
println(s)
break
end
i+=1
x0=x
end
s=@sprintf("Fixed-point failed after %d iteration.\n",i)
println(s)

輸出畫面

0    0.0000000
1    0.5773503
2    0.7705652
3    0.8487220
4    0.8825453
5    0.8975975
6    0.9043784
7    0.9074499
8    0.9088446
9    0.9094786
10   0.9097669
11   0.9098981
The Root= 0.9099578  x-x0= 0.0000597
Fixed-point failed after 12 iteration.

0    0.0000000
1   -0.5773503
2   -0.4325829
3   -0.4650559
4   -0.4575660
5   -0.4592828
6   -0.4588887
The Root=-0.4589791  x-x0= 0.0000904
Fixed-point failed after 7 iteration.

Julia語言例題 EX2-9 定點回路法求非線性方程式 f(x)=1/5^x - x = 0

Julia語言例題 EX2-9 定點回路法求非線性方程式 f(x)=1/5^x - x = 0

F(X)= =1/5^x - x 改寫成
g(x)= 1/ (5^x) 與 g(x)= x
f(0.4) * f(0.5) < 0 所以取 x0=0.45

using Printf
#=======================================================
* ex2-9.jl is used for solving nonlinear equation
* based on Fixed-Point Algorithm g(x)=x with initial
* approximation P0.
=========================================================#
MAX=50
TOL=0.0001

function g(x0::Float64)
return (1/(5^x))
end

i=1
x0=0.45
x=0.0
while(i<=MAX)
x=g(x0)
s=@sprintf("%-2d %10.7lf",i-1,x0)
println(s)
if(abs(x-x0) < TOL)
s=@sprintf("The Root=%10.7lf x-x0=%10.7lf",x,abs(x-x0))
println(s)
break
end
i+=1
x0=x
end
s=@sprintf("Fixed-point failed after %d iteration.\n",i)
println(s)

輸出結果

0    0.4500000
1    1.0000000
2    0.2000000
3    0.7247797
4    0.3114589
5    0.6057586
6    0.3772185
7    0.5449236
8    0.4160205
9    0.5119342
10   0.4387058
11   0.4935803
12   0.4518582
13   0.4832420
14   0.4594395
15   0.4773815
16   0.4637935
17   0.4740479
18   0.4662885
19   0.4721482
20   0.4677164
21   0.4710644
22   0.4685329
23   0.4704457
24   0.4689997
25   0.4700925
26   0.4692664
27   0.4698907
28   0.4694188
29   0.4697755
30   0.4695059
31   0.4697096
32   0.4695556
33   0.4696720
The Root= 0.4695841  x-x0= 0.0000880
Fixed-point failed after 34 iteration.

Julia語言範例2-5 非線性方程式 f(x)= (4x-7) /(x-1)的根

Julia語言範例2-5 非線性方程式 f(x)= (4x-7) /(x-1)
利用牛頓(Newton-Raphson)方法求f(x)的根取誤差=0.001
x0=1.5 ,1.625 , 1.875 , 1.95 , 3.0 時的現象 (會出現overflow)

using Printf
#=========================================================
Julia has no do-while construct. Here is one of several
ways to implement do-while behavior.

julia> i = 0
0

julia> while true
println(i)
i += 1
i % 6 == 0 && break
end
============================================================#

MAX=100
TOL=0.001

function f(x0::Float64)
return ((4*x0-7)/(x0-2))
end

function ff(x0::Float64)
return (-1/(x0-2)^2)
end

i=1
x0=1.5
while(i<=MAX)
x=x0-f(x0)/ff(x0)
s=@sprintf("%2d %10.7lf",i,x0)
println(s)
if(abs(x-x0)<TOL)
s=@sprintf("Root=%10.7lf x-x0=%10.7lf\n",x,abs(x-x0))
println(s)
break
end
i+=1
x0=x
end
s=@sprintf("Newton-Raphson Method failed after %2d iterations!!!\n",i)
println(s)

i=1
x0=1.625
while(i<=MAX)
x=x0-f(x0)/ff(x0)
s=@sprintf("%2d %10.7lf",i,x0)
println(s)
if(abs(x-x0)<TOL)
s=@sprintf("Root=%10.7lf x-x0=%10.7lf\n",x,abs(x-x0))
println(s)
break
end
i+=1
x0=x
end
s=@sprintf("Newton-Raphson Method failed after %2d iterations!!!\n",i)
println(s)

i=1
x0=1.875
while(i<=MAX)
x=x0-f(x0)/ff(x0)
s=@sprintf("%2d %10.7lf",i,x0)
println(s)
if(abs(x-x0)<TOL)
s=@sprintf("Root=%10.7lf x-x0=%10.7lf\n",x,abs(x-x0))
println(s)
break
end
i+=1
x0=x
end
s=@sprintf("Newton-Raphson Method failed after %2d iterations!!!\n",i)
println(s)

i=1
x0=1.95
while(i<=MAX)
x=x0-f(x0)/ff(x0)
s=@sprintf("%2d %10.7lf",i,x0)
println(s)
if(abs(x-x0)<TOL)
s=@sprintf("Root=%10.7lf x-x0=%10.7lf\n",x,abs(x-x0))
println(s)
break
end
i+=1
x0=x
end
s=@sprintf("Newton-Raphson Method failed after %2d iterations!!!\n",i)
println(s)

i=1
x0=3.0
while(i<=MAX)
x=x0-f(x0)/ff(x0)
s=@sprintf("%2d %10.7lf",i,x0)
println(s)
if(abs(x-x0)<TOL)
s=@sprintf("Root=%10.7lf x-x0=%10.7lf\n",x,abs(x-x0))
println(s)
break
end
i+=1
x0=x
end
s=@sprintf("Newton-Raphson Method failed after %2d iterations!!!\n",i)
println(s)

輸出畫面

x0= 1.5  Newton-Raphson Method failed after 101 iterations!!!

x0=1.625 Root= 1.7500038 x-x0= 0.0009727
x0=1.875 Root= 1.7500038 x-x0= 0.0009727

x0=1.95  Root= 1.7500002 x-x0= 0.0001979

x0=3.0   Newton-Raphson Method failed after 101 iterations!!!

======================================================

 1    1.5000000
 2    2.0000000
 3          NaN
 4          NaN
 5          NaN
 6          NaN
 7          NaN
 8          NaN
 9          NaN
10          NaN
11          NaN
12          NaN
13          NaN
14          NaN
15          NaN
16          NaN
17          NaN
18          NaN
19          NaN
20          NaN
21          NaN
22          NaN
23          NaN
24          NaN
25          NaN
26          NaN
27          NaN
28          NaN
29          NaN
30          NaN
31          NaN
32          NaN
33          NaN
34          NaN
35          NaN
36          NaN
37          NaN
38          NaN
39          NaN
40          NaN
41          NaN
42          NaN
43          NaN
44          NaN
45          NaN
46          NaN
47          NaN
48          NaN
49          NaN
50          NaN
51          NaN
52          NaN
53          NaN
54          NaN
55          NaN
56          NaN
57          NaN
58          NaN
59          NaN
60          NaN
61          NaN
62          NaN
63          NaN
64          NaN
65          NaN
66          NaN
67          NaN
68          NaN
69          NaN
70          NaN
71          NaN
72          NaN
73          NaN
74          NaN
75          NaN
76          NaN
77          NaN
78          NaN
79          NaN
80          NaN
81          NaN
82          NaN
83          NaN
84          NaN
85          NaN
86          NaN
87          NaN
88          NaN
89          NaN
90          NaN
91          NaN
92          NaN
93          NaN
94          NaN
95          NaN
96          NaN
97          NaN
98          NaN
99          NaN
100          NaN
Newton-Raphson Method failed after 101 iterations!!!

 1    1.6250000
 2    1.8125000
 3    1.7656250
 4    1.7509766
Root= 1.7500038 x-x0= 0.0009727

Newton-Raphson Method failed after  4 iterations!!!

 1    1.8750000
 2    1.8125000
 3    1.7656250
 4    1.7509766
Root= 1.7500038 x-x0= 0.0009727

Newton-Raphson Method failed after  4 iterations!!!

 1    1.9500000
 2    1.9100000
 3    1.8524000
 4    1.7919430
 5    1.7570369
 6    1.7501981
Root= 1.7500002 x-x0= 0.0001979

Newton-Raphson Method failed after  6 iterations!!!

 1    3.0000000
 2    8.0000000
 3   158.0000000
 4   97658.0000000
 5   38146972658.0000229
 6   5820766091346748375040.0000000
 7   135525271560688433474159118908309231747727360.0000000
 8   73468396926393384679514105682249657134065937588368498667350306457576200536943192421957632.0000000
 9   21590421387736356954845175477318195813934831457681903459410459099972442985296985166524414010377109304944487332172036675687278085164774782269901823598293450947135537747747397435392.0000000
10          Inf
11          NaN
12          NaN
13          NaN
14          NaN
15          NaN
16          NaN
17          NaN
18          NaN
19          NaN
20          NaN
21          NaN
22          NaN
23          NaN
24          NaN
25          NaN
26          NaN
27          NaN
28          NaN
29          NaN
30          NaN
31          NaN
32          NaN
33          NaN
34          NaN
35          NaN
36          NaN
37          NaN
38          NaN
39          NaN
40          NaN
41          NaN
42          NaN
43          NaN
44          NaN
45          NaN
46          NaN
47          NaN
48          NaN
49          NaN
50          NaN
51          NaN
52          NaN
53          NaN
54          NaN
55          NaN
56          NaN
57          NaN
58          NaN
59          NaN
60          NaN
61          NaN
62          NaN
63          NaN
64          NaN
65          NaN
66          NaN
67          NaN
68          NaN
69          NaN
70          NaN
71          NaN
72          NaN
73          NaN
74          NaN
75          NaN
76          NaN
77          NaN
78          NaN
79          NaN
80          NaN
81          NaN
82          NaN
83          NaN
84          NaN
85          NaN
86          NaN
87          NaN
88          NaN
89          NaN
90          NaN
91          NaN
92          NaN
93          NaN
94          NaN
95          NaN
96          NaN
97          NaN
98          NaN
99          NaN
100          NaN
Newton-Raphson Method failed after 101 iterations!!!

julia語言牛頓法(Newton - Raphson)以圖形概念主要是不停在取切線斜率。

julia語言牛頓法以圖形概念主要是不停在取切線斜率。

虛擬碼
Algorithm NewtonRoot

E0 : 初始化最小誤差 EPS，初始點 xo

E1 : x = x0

E2 : x0 = x - f(x) / f'(x)

E3 : if abs(x-x0) < eps，演算法結束，傳回 x

E4 : goto E1

End Algorithm

#===========================================================
> func(-1.781503813804761e+000) = +0.000000000000000e+000
> func(+2.313837835588586e+000) = +0.000000000000000e+000
> func(+4.947665978216175e+000) = +3.552713678800501e-015
> func(-1.781503813804761e+000) = +0.000000000000000e+000

Julia has no do-while construct. Here is one of several ways to implement do-while behavior.

julia> i = 0
0

julia> while true
println(i)
i += 1
i % 6 == 0 && break
end

===========================================================#

using Printf

#[ -2.00 , -1.00 ] , [ 2.00 , 3.00 ] , [ +4.00 , +5.00 ]

function func(x::Float64)
x2=x*x
x3=x2*x;
return (x3 - 5.48*x2 - 1.4883*x + 20.394828)
end

# funcd(x) = func'(x)
function funcd(x::Float64)
x2=x*x
return (3*x2-10.96*x-1.4883)
end

#-------------------------------------------------------

function NewtonRoot(x0::Float64,eps::Float64) #/* 初點容許誤差*/

x=x0;
while true
x0=x;
x = x0 - func(x0) / funcd(x0)
(abs(x-x0)<eps) && break
end

return x
end

eps=1.0E-9
max_iterator=100

x0 = -2.0
y0 = NewtonRoot(x0, eps)
s=@sprintf("\n> func(%+.15e) = %+.15e", y0, func(y0))
println(s)

x0 = 2.0
y0 = NewtonRoot(x0, eps)
s=@sprintf("\n> func(%+.15e) = %+.15e", y0, func(y0))
println(s)

x0 = 4.0
y0 = NewtonRoot(x0, eps)
s=@sprintf("\n> func(%+.15e) = %+.15e", y0, func(y0))
println(s)

x0 = -10.0 #// test
y0 = NewtonRoot(x0, eps)
s=@sprintf("\n> func(%+.15e) = %+.15e", y0, func(y0))
println(s)

輸出畫面

> func(-1.781503813804761e+00) = +0.000000000000000e+00

> func(+2.313837835588586e+00) = +0.000000000000000e+00

> func(+4.947665978216175e+00) = +3.552713678800501e-15

> func(-1.781503813804761e+00) = +0.000000000000000e+00

Julia語言例題2-6 已知方程式 e^x + x^-2 + 2 cosx -6 利用正割法找出f(x)=0的根

Julia語言例題2-6 已知方程式 e^x + x^-2 + 2 cosx -6 利用正割法找出f(x)=0的根

虛擬碼
Algorithm SecantRoot

E0 : 初始化最小誤差 EPS，初始點 xo, x1

E1 : 算 y0 = f(x0), y1 = f(x1)

E2 : x2 = x1 - y1*(x1-x0) / (y1-y0)， delta = x2-x1

E3 : x0 = x1, x1 = x2, y0=y1, y1=f(x1)

E4 : if abs(delta) < eps 或 abs(y1-y0) < eps，演算法結束，傳回 x2

E5 : goto E2

End Algorithm

using Printf
#========================================================================
/*----------------------------------------------------------------*\
|
| E0 : 初始化最小誤差EPS，初始點x0, x1
| E1 : y0 = f(x0) , y1 = f(x1)
| E2 : x2 = x1 - y1 * (x1 - x0) / (y1 - y0), delta = x2 - x1
| E3 : x0 = x1, x1 = x2, y0 = y1, y1=f(x1)
| E4 : if abs(delta) or abs(y1-y0), 演算法結束, 傳回x2
| E5 : goto E2
|
\*----------------------------------------------------------------*/
========================================================================#

MAX=50
TOL=0.00001

function fx(x::Float64)
return (exp(x)+1/(2^x)+2*cos(x)-6)
end

i=2
x0=1.8
x1=2.0
q0=fx(x0);
q1=fx(x1);
s=@sprintf("i xi f(x)")
println(s)
s=@sprintf("%-2d %10.6lf %10.6lf",0,x0,q0)
println(s)
s=@sprintf("%-2d %10.6lf %10.6lf",1,x1,q1)
println(s)

while(i<=MAX)
x=x1-q1*(x1-x0)/(q1-q0);
s=@sprintf("%-2d %10.6lf %10.6lf",i,x,fx(x))
println(s)
if(abs(x-x1) < TOL)
s=@sprintf("The Root=%10.6lf f(%10.6lf)=%10.6lf",x,x,fx(x))
println(s)
break
else
i+=1
x0=x1
q0=q1
x1=x
q1=fx(x)
end
end

if(i>MAX)
s=@printf("Secant Method faileds!!!\n")
println(s)
end

輸出畫面

i       xi           f(x)
0      1.800000    -0.117582
1      2.000000     0.806762
2      1.825441    -0.016116
3      1.828860    -0.002147
4      1.829385     0.000007
5      1.829384    -0.000000
The Root=  1.829384    f(  1.829384)= -0.000000

非線性方程式求解 - 目錄

非線性方程式求解 - 目錄

源自於 http://edisonx.pixnet.net/blog/post/83349249

[C語言數值分析] 非線性方程式求解 - 勘根定理

[C語言數值分析] 非線性方程式求解 - 二分法 BiSect

[C語言數值分析] 非線性方程式求解 - 牛頓法 Newton-Raphson

[C語言數值分析] 非線性方程式求解 - 割線法 Secant

[C語言數值分析] 非線性方程式求解 - 假位法 False Position

[C語言數值分析] 非線性方程式求解 - 定點法 Fix Point

[C語言數值分析] 非線性方程式求解 - 艾肯疊代法 Aitken

[C語言數值分析] 非線性方程式求解 - 穆勒法 Muller

[C語言數值分析] 非線性方程式求解 - Force Find a root

非線性方程式求解 - 目錄

Julia語言範例2-3 使用Newton-Raphson理則解 x^3 + 4x^2 -10 =0 非線性方程式解

Julia語言範例2-3 使用Newton-Raphson理則
解 x^3 + 4x^2 -10 =0 非線性方程式解
使用Newton-Raphson理則找出位於(1.0 ,2.0)之間的根
取誤差=0.001

虛擬碼

Algorithm NewtonRoot

E0 : 初始化最小誤差 EPS，初始點 xo

E1 : x = x0

E2 : x0 = x - f(x) / f'(x)

E3 : if abs(x-x0) < eps，演算法結束，傳回 x

E4 : goto E1

End Algorithm

==========================================
using Printf

MAX=100
TOL=0.001
function fx(x::Float64)
return (x^3 + 4*x^2 -10)
end

function ffx(x::Float64)
return (3*x^2 + 8x)
end

i=1
x0=1.5
println(" i Root of f(x)")
while(i<=MAX)
x=x0-fx(x0)/ffx(x0);
s=@sprintf("%2d %10.7lf",i,x0);
println(s)
if(abs(x-x0) <TOL)
s=@sprintf("Root=%10.7lf x-x0=%10.7lf",x,abs(x-x0));
println(s)
break
end
i+=1
x0=x
end
println()
s=@sprintf("Newton-Raphson Method failed after %2d iterations!!!\n",i)
println(s)

輸出畫面

 i     Root of f(x)
 1    1.5000000
 2    1.3733333
 3    1.3652620
Root= 1.3652300 x-x0= 0.0000320

Newton-Raphson Method failed after  3 iterations!!!

Julia語言範例2-4 使用Newton-Raphson理則解 cos(x) - x =0 非線性方程式解

Julia語言例題2-4 使用Newton-Raphson理則解 cos(x) - x =0 非線性方程式解

using Printf

MAX=100
TOL=0.001
function fx(x::Float64)
return (cos(x)-x)
end

function ffx(x::Float64) # fx =(cos(x)-x) 的微分值
return (-sin(x)-1)
end

i=1
x0=1.0
println(" i Root of f(x)")
while(i<=MAX)
x=x0-fx(x0)/ffx(x0);
s=@sprintf("%2d %10.7lf",i,x0);
println(s)
if(abs(x-x0) <TOL)
s=@sprintf("Root=%10.7lf x-x0=%10.7lf",x,abs(x-x0));
println(s)
break
end
i+=1
x0=x
end
println()
s=@sprintf("Newton-Raphson Method failed after %2d iterations!!!\n",i)
println(s)

輸出畫面

 i     Root of f(x)
 1    1.0000000
 2    0.7503639
 3    0.7391129
Root= 0.7390851 x-x0= 0.0000278

Newton-Raphson Method failed after  3 iterations!!!

Julia語言 EX2-2 f(x) = (0.0000115pow(x,2)+0.000014pow(x,1.5)-0.01962) 利用二分法求解---(a) 測試根的範圍

Julia語言 EX2-2 f(x) = (0.0000115*pow(x,2)+0.000014*pow(x,1.5)-0.01962) 利用二分法求解---(a) 測試根的範圍

using Printf

function f(x::Float64)
return (0.0000115*(x^2) + 0.000014*(x^1.5)- 0.01962)
end

#define MAX 50 /* maximum iterations */
#define TOL 0.001 /* tolerance */
MAX=150
TOL=0.001
i=1
a=0.0
b=50.0
x=a
s=@sprintf(" i\t\tx\t\t\tf(x)")
println(s)
while(x<=b)
s=@sprintf("%2d\t\t%9.6lf\t\t%9.6lf", i,x,f(x))
println(s)
i+=1
x+=0.1
end

a=0.0 , b=50.0的範圍
輸出畫面

 i      x      f(x)
 1   0.000000  -0.019620
 2   0.100000  -0.019619
................................................

140  13.900000  -0.016673
141  14.000000  -0.016633
142  14.100000  -0.016592
143  14.200000  -0.016552
144  14.300000  -0.016511
145  14.400000  -0.016470
146  14.500000  -0.016429
................................................
357  35.600000  -0.002072
358  35.700000  -0.001977
359  35.800000  -0.001882
360  35.900000  -0.001787
361  36.000000  -0.001692
362  36.100000  -0.001596
363  36.200000  -0.001501
364  36.300000  -0.001405
365  36.400000  -0.001308
366  36.500000  -0.001212
367  36.600000  -0.001115
368  36.700000  -0.001018
369  36.800000  -0.000921
370  36.900000  -0.000823
371  37.000000  -0.000726
372  37.100000  -0.000628
373  37.200000  -0.000529
374  37.300000  -0.000431
375  37.400000  -0.000332
376  37.500000  -0.000233
377  37.600000  -0.000134
378  37.700000  -0.000034
379  37.800000   0.000065
380  37.900000   0.000165
381  38.000000   0.000265
382  38.100000   0.000366
383  38.200000   0.000467
384  38.300000   0.000568
385  38.400000   0.000669
386  38.500000   0.000770
387  38.600000   0.000872

.................................................

493  49.200000   0.013049

494  49.300000   0.013177
495  49.400000   0.013305
496  49.500000   0.013434
497  49.600000   0.013562
498  49.700000   0.013691
499  49.800000   0.013821

500 49.900000 0.013950

2019年2月24日星期日

Julia語言 EX2-2 f(x)=(0.0000115pow(x,2)+0.000014pow(x,1.5)-0.01962) 二分法求解

Julia語言 EX2-2 f(x)=(0.0000115*pow(x,2)+0.000014*pow(x,1.5)-0.01962) 二分法求解誤差值 TOL=0.001

using Printf

function f(x::Float64)
return (0.0000115*(x^2) + 0.000014*(x^1.5)- 0.01962)
end

#define MAX 50 /* maximum iterations */
#define TOL 0.001 /* tolerance */
MAX=50
TOL=0.001
i=1
a=37.0
b=38.0
s=@sprintf(" i\t\ta(i)\t\t\tb(i)\t\t\tp(i)\t\t\tf(p(i))\t\t\th(i)\n")
println(s)
while(i<=MAX)
p=(a+b)/2.0
h=abs((b-a)/2.0)
s=@sprintf("%2d\t\t%9.6lf\t\t%9.6lf\t\t%9.6lf\t\t%9.6lf\t\t%9.6lf", i,a,b,p,f(p),h)
println(s)
if(abs(f(p))==0 || h<TOL)
s=@sprintf("%2d iterations!!!\n",i);
println(s)
S=@sprintf(" The Root=%9.6lf f(%9.6lf)=%9.6lf",p,p,f(p));
println(s)
break
else
if(f(a)*f(p) >0)
a=p
i+=1
else(f(b)*f(p) >0)
b=p
i+=1
end
end

end
s=@sprintf("Bisection Method failed after %d iterations!!!\n",i);
println(s)

輸出畫面

 i  a(i)   b(i)   p(i)   f(p(i))   h(i)

 1  37.000000  38.000000  37.500000  -0.000233   0.500000
 2  37.500000  38.000000  37.750000   0.000015   0.250000
 3  37.500000  37.750000  37.625000  -0.000109   0.125000
 4  37.625000  37.750000  37.687500  -0.000047   0.062500
 5  37.687500  37.750000  37.718750  -0.000016   0.031250
 6  37.718750  37.750000  37.734375  -0.000000   0.015625
 7  37.734375  37.750000  37.742188   0.000008   0.007813
 8  37.734375  37.742188  37.738281   0.000004   0.003906
 9  37.734375  37.738281  37.736328   0.000002   0.001953
10  37.734375  37.736328  37.735352   0.000001   0.000977
10 iterations!!!

10 iterations!!!

Bisection Method failed after 10 iterations!!!

Julia語言 EX2-2 f(x) = (0.0000115pow(x,2)+0.000014pow(x,1.5)-0.01962) 利用二分法求解

Julia語言 EX2-2 f(x) = (0.0000115*pow(x,2)+0.000014*pow(x,1.5)-0.01962) 利用二分法求解

using Printf

#==============================================================
Arithmetic Operators
The following arithmetic operators are supported on all primitive numeric types:

Expression Name Description
+x unary plus the identity operation
-x unary minus maps values to their additive inverses
x + y binary plus performs addition
x - y binary minus performs subtraction
x * y times performs multiplication
x / y divide performs division
x ÷ y integer divide x / y, truncated to an integer
x \ y inverse divide equivalent to y / x
x ^ y power raises x to the yth power
x % y remainder equivalent to rem(x,y)
==============================================================#

function fx(x::Float64) # 原非線性方程式
#( (1.15*10^(-5)*(x^2)) + (1.4*10^(-5)*x) - (1.962*10^(-2)) )
y= (0.0000115*(x^2)+(0.000014*(x^1.5)) - 0.1962)
return ( y )
end

a=100.0
b=200.0
x=a
println(" i\t x\t f(x)");
i=1
while (x<=b)
s=@sprintf("%2d\t%5.2f\t%10.5f",i,x,fx(x))
println(s)
x=x+1.0
i+=1
end

a=124.0
b=125.0
x=a
println(" i\t x\t f(x)");
i=1
while (x<=b)
s=@sprintf("%2d\t%5.2f\t%10.5f",i,x,fx(x))
println(s)
x=x+0.001
i+=1
end

i   x    f(x)
 1 100.00   -0.06720
 2 101.00   -0.06468
 3 102.00   -0.06213
 4 103.00   -0.05956
 5 104.00   -0.05697
 6 105.00   -0.05435
 7 106.00   -0.05171
 8 107.00   -0.04904
 9 108.00   -0.04635
10 109.00   -0.04364
11 110.00   -0.04090
12 111.00   -0.03814
13 112.00   -0.03535
14 113.00   -0.03254
15 114.00   -0.02971
16 115.00   -0.02685
17 116.00   -0.02396
18 117.00   -0.02106
19 118.00   -0.01813
20 119.00   -0.01517
21 120.00   -0.01220
22 121.00   -0.00919
23 122.00   -0.00617
24 123.00   -0.00312
25 124.00   -0.00004
26 125.00    0.00305
27 126.00    0.00617
28 127.00    0.00932
29 128.00    0.01249
30 129.00    0.01568
31 130.00    0.01890
32 131.00    0.02214
33 132.00    0.02541
34 133.00    0.02870
35 134.00    0.03201
36 135.00    0.03535
37 136.00    0.03871
38 137.00    0.04209
39 138.00    0.04550
40 139.00    0.04893
41 140.00    0.05239
42 141.00    0.05587
43 142.00    0.05938
44 143.00    0.06290
45 144.00    0.06646
46 145.00    0.07003
47 146.00    0.07363
48 147.00    0.07726
49 148.00    0.08090
50 149.00    0.08457
51 150.00    0.08827
52 151.00    0.09199
...................
94 193.00    0.26970
95 194.00    0.27444
96 195.00    0.27921
97 196.00    0.28400
98 197.00    0.28881
99 198.00    0.29365
100 199.00    0.29851
101 200.00    0.30340

 i   x    f(x)
 1 124.00   -0.00004
 2 124.00   -0.00004
 3 124.00   -0.00004
 4 124.00   -0.00004
 5 124.00   -0.00003
 6 124.01   -0.00003
 7 124.01   -0.00003
 8 124.01   -0.00002
 9 124.01   -0.00002
10 124.01   -0.00002
11 124.01   -0.00001
12 124.01   -0.00001
13 124.01   -0.00001
14 124.01   -0.00000
15 124.01   -0.00000
16 124.02    0.00000
17 124.02    0.00000
18 124.02    0.00001
19 124.02    0.00001
20 124.02    0.00001
21 124.02    0.00002
22 124.02    0.00002
23 124.02    0.00002
24 124.02    0.00003
25 124.02    0.00003
26 124.03    0.00003
27 124.03    0.00004
28 124.03    0.00004
29 124.03    0.00004
30 124.03    0.00004
31 124.03    0.00005
32 124.03    0.00005
33 124.03    0.00005
34 124.03    0.00006
35 124.03    0.00006
36 124.04    0.00006
37 124.04    0.00007
38 124.04    0.00007
39 124.04    0.00007
40 124.04    0.00008
41 124.04    0.00008
42 124.04    0.00008
43 124.04    0.00008
44 124.04    0.00009
45 124.04    0.00009
46 124.05    0.00009
47 124.05    0.00010
48 124.05    0.00010
49 124.05    0.00010
50 124.05    0.00011
51 124.05    0.00011
52 124.05    0.00011
53 124.05    0.00012
54 124.05    0.00012
55 124.05    0.00012
56 124.06    0.00013
57 124.06    0.00013
58 124.06    0.00013
59 124.06    0.00013

..........................

994 124.99    0.00303
995 124.99    0.00303
996 125.00    0.00304
997 125.00    0.00304
998 125.00    0.00304
999 125.00    0.00305
1000 125.00    0.00305

Julia語言 EX2-1利用2分法求 x^3-4x+2 =0 在0.0到2.0之間的解

Julia語言 EX2-1利用2分法求 x^3-4x+2 =0 在0.0到2.0之間的解
================================
x^3-4x+2 =0 在0.0到2.0之間的解
================================

虛擬碼

Algorithm BiSect

E0 : 初始化最小誤差 EPS，初始化上界 a，下界 b。

E1 : c = (a+b) / 2

E2 : if abs ( f(c) ) <= EPS ，演算法結束，傳回 c 值。

E3 : if ( f(c) * f(a) <= 0) b = c;
else a = c;

E4 : goto E1

End Algorithm

程式

using Printf

#==============================================================
Arithmetic Operators
The following arithmetic operators are supported on all primitive numeric types:

Expression Name Description
+x unary plus the identity operation
-x unary minus maps values to their additive inverses
x + y binary plus performs addition
x - y binary minus performs subtraction
x * y times performs multiplication
x / y divide performs division
x ÷ y integer divide x / y, truncated to an integer
x \ y inverse divide equivalent to y / x
x ^ y power raises x to the yth power
x % y remainder equivalent to rem(x,y)
==============================================================#

function fx(x::Float64) # 原非線性方程式
return (x^3 -4.0*x +2.0)
end

a=0.0
b=2.0
x=a
println(" i\t x\t f(x)");
i=1
while (x<=b)
s=@sprintf("%2d\t%5.2f\t%10.5f",i,x,fx(x))
println(s)
x=x+0.1
i+=1
end

輸出畫面

 i   x    f(x)
 1  0.00    2.00000
 2  0.10    1.60100
 3  0.20    1.20800
 4  0.30    0.82700
 5  0.40    0.46400
 6  0.50    0.12500
 7  0.60   -0.18400
 8  0.70   -0.45700
 9  0.80   -0.68800
10  0.90   -0.87100
11  1.00   -1.00000
12  1.10   -1.06900
13  1.20   -1.07200
14  1.30   -1.00300
15  1.40   -0.85600
16  1.50   -0.62500
17  1.60   -0.30400
18  1.70    0.11300
19  1.80    0.63200
20  1.90    1.25900

Julia語言習題1-7 已知下面4點座標求Hermite內插法的插除表並計算Hn(x) = [1.5 , 2.5 , 3.5 , 4.5]之值

Julia語言習題1-7 已知下面4點座標求Hermite內插法的插除表並計算Hn(x) = [1.5 , 2.5 , 3.5 , 4.5]之值

已知下面4點座標

x f(x) f'(x)=ff(x)
=================
1.0 0.000 1.0000
2.0 0.693 0.5000
3.0 1.099 0.3333

4.0 1.386 0.2500
==================
若已知 f(x) = ln(x)
求Hermite內插法的插除表並
計算Hn(x) = [1.5 , 2.5 , 3.5 , 4.5]之值

using Printf

#======================================================

function [H]=hermite(X,x,f,fd);
m=length(x);
for i=1:m
z(2*i-1)=x(i);
Q(2*i-1,1)=f(i);
z(2*i)=x(i);
Q(2*i,1)=f(i);
Q(2*i,2)=fd(i);
if i~=1
Q(2*i-1,2)=(Q(2*i-1,1)-Q(2*i-2,1))/(z(2*i-1)-z(2*i-2));
end;
end;
for i=3:2*m
for j=3:i
Q(i,j)=(Q(i,j-1)-Q(i-1,j-1))/(z(i)-z(i-j+1));
end
end

p=1;
H=Q(1,1);
for i=2:2*m
p=p.*(X-z(i-1));
H=H+p.*Q(i,i);
end;

======================================================#
println("Hermite 差除表")
x=[ 1.0 ,2.0 ,3.0 ,4.0 ]
f=[0.0 , 0.693 , 1.099 , 1.386]
ff=[1.00 , 0.5 , 0.3333 , 0.25 ]

xx=[1.5 , 2.5 , 3.5 , 4.5]
n=length(x)
z=[0.0 for i=1:2*n]

# q 需為 f 的 2n * 2n 倍
q= [[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ],
[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ],
[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ],
[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ],
[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ],
[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ],
[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ],
[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.0 ]]

println("i\tz(i)\tf(i)\t\tf(i-1,i)\tf(i-2,i-1,i)...\n");

for i=1:n
z[2*i-1]=x[i]
z[2*i]=x[i]

q[2*i-1][1]=f[i]
q[2*i][1]=f[i]
q[2*i][2]=ff[i]

if (i!=1)
q[2*i-1][2]= (q[2*i-1][1]- q[2*i-2][1]) / (z[2*i-1]-z[2*i-2])
end
#println(i,z,q)
end

for i=1:n
if(i==1)
s=@sprintf("%d\t%0.4f\t%0.4f\n",i,z[i],q[i][1])
print(s)
else (i==2)
s=@sprintf("%d\t%0.4f\t%0.4f\t\t%0.4f\n",i,z[i],q[i][1],q[i][2])
print(s)
end
end

for i=3:2*n
s=@sprintf("%d\t%0.4f\t%0.4f\t\t%0.4f",i,z[i],q[i][1],q[i][2])
print(s)
for j=3:i
q[i][j]=(q[i][j-1]-q[i-1][j-1])/(z[i]-z[i-j+1])
s=@sprintf("\t%0.4f",q[i][j])
print(s)
end
println()
end

#println(q)

for m=1:length(xx)
xa=xx[m]
println()
println("============================================")
println("i=",m," xa=",xa)

p=1
H=q[1][1]
for i=2:2*n
p *= (xa-z[i-1])
H = H+p*q[i][i]
end
println()

println("Xa=",xa," Hn(xa)=" ,H, "---Hermite方法")
H1=(log(xa))
println("真實解 f(x)= ln(x)")
println("Xa=",xa," f(xa)=" ,H1 , "---真實值")
println("Hn(x)與fn(x)的誤差" ,abs(H1-H))
end

輸出畫面

Hermite 差除表
i z(i) f(i)  f(i-1,i) f(i-2,i-1,i)...

1 1.0000 0.0000
2 1.0000 0.0000  1.0000
3 2.0000 0.6930  0.6930
4 2.0000 0.6930  0.5000
3 2.0000 0.6930  0.6930 -0.3070
4 2.0000 0.6930  0.5000 -0.1930 0.1140
5 3.0000 1.0990  0.4060 -0.0940 0.0495 -0.0323
6 3.0000 1.0990  0.3333 -0.0727 0.0213 -0.0141 0.0091
7 4.0000 1.3860  0.2870 -0.0463 0.0132 -0.0040 0.0034 -0.0019
8 4.0000 1.3860  0.2500 -0.0370 0.0093 -0.0019 0.0011 -0.0008 0.0004

============================================
i=1  xa=1.5

Xa=1.5 Hn(xa)=0.4057314453125---Hermite方法
真實解 f(x)= ln(x)
Xa=1.5 f(xa)=0.4054651081081644---真實值
Hn(x)與fn(x)的誤差0.00026633720433560937

============================================
i=2  xa=2.5

Xa=2.5 Hn(xa)=0.9164583984375---Hermite方法
真實解 f(x)= ln(x)
Xa=2.5 f(xa)=0.9162907318741551---真實值
Hn(x)與fn(x)的誤差0.0001676665633448815

============================================
i=3  xa=3.5

Xa=3.5 Hn(xa)=1.2529150390625001---Hermite方法
真實解 f(x)= ln(x)
Xa=3.5 f(xa)=1.252762968495368---真實值
Hn(x)與fn(x)的誤差0.00015207056713206768

============================================
i=4  xa=4.5

Xa=4.5 Hn(xa)=1.5076044921874998---Hermite方法
真實解 f(x)= ln(x)
Xa=4.5 f(xa)=1.5040773967762742---真實值
Hn(x)與fn(x)的誤差0.0035270954112256447

Julia語言例題1-12 已知3點座標及導數 ,建立 Hermite 向後差除表

Julia語言例題1-12 已知3點座標及導數 ,建立 Hermite 向後差除表

例題1-12 已知3點座標及導數如下

=================================
i xi f(xi) f'(xi)
0 0.0 1.0000 1.0000
1 1.0 2.7183 2.7183
2 2.0 2.7183 2.7183
==================================
建立 Hermite 向後差除表
Hn(1.5) 之值 , H'n(1.5) 之值

程式
using Printf

#======================================================

function [H]=hermite(X,x,f,fd);
m=length(x);
for i=1:m
z(2*i-1)=x(i);
Q(2*i-1,1)=f(i);
z(2*i)=x(i);
Q(2*i,1)=f(i);
Q(2*i,2)=fd(i);
if i~=1
Q(2*i-1,2)=(Q(2*i-1,1)-Q(2*i-2,1))/(z(2*i-1)-z(2*i-2));
end;
end;
for i=3:2*m
for j=3:i
Q(i,j)=(Q(i,j-1)-Q(i-1,j-1))/(z(i)-z(i-j+1));
end
end

p=1;
H=Q(1,1);
for i=2:2*m
p=p.*(X-z(i-1));
H=H+p.*Q(i,i);
end;

======================================================#
println("Hermite 差除表")
x=[0.0 , 1.0 ,2.0]
f=[1.00 , 2.7183 , 7.3891]
ff=[1.00 , 2.7183 , 7.3891]

n=length(x)
z=[0.0 for i=1:2*n]

# q 需為 f 的 2n * 2n 倍
q= [[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ],
[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ],
[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ],
[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ],
[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ],
[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ]]

println("i\tz(i)\tf(i)\t\tf(i-1,i)\tf(i-2,i-1,i)...\n");

for i=1:n
z[2*i-1]=x[i]
z[2*i]=x[i]

q[2*i-1][1]=f[i]
q[2*i][1]=f[i]
q[2*i][2]=ff[i]

if (i!=1)
q[2*i-1][2]= (q[2*i-1][1]- q[2*i-2][1]) / (z[2*i-1]-z[2*i-2])
end
#println(i,z,q)
end

for i=1:n
if(i==1)
s=@sprintf("%d\t%0.4f\t%0.4f\n",i,z[i],q[i][1])
print(s)
else (i==2)
s=@sprintf("%d\t%0.4f\t%0.4f\t\t%0.4f\n",i,z[i],q[i][1],q[i][2])
print(s)
end
end

for i=3:2*n
s=@sprintf("%d\t%0.4f\t%0.4f\t\t%0.4f",i,z[i],q[i][1],q[i][2])
print(s)
for j=3:i
q[i][j]=(q[i][j-1]-q[i-1][j-1])/(z[i]-z[i-j+1])
s=@sprintf("\t%0.4f",q[i][j])
print(s)
end
println()
end

#println(q)

xa=1.5
p=1
H=q[1][1]

for i=2:2*n
p *= (xa-z[i-1])
H = H+p*q[i][i]
end
println()

println("Xa=",xa," Hn(xa)=" ,H, "---Hermite方法")
H1=(exp(xa))
println("真實解 f(x)=exp(x)")
println("Xa=",xa," f(xa)=" ,H1 , "---真實值")
println("Hn(x)與fn(x)的誤差" ,abs(H1-H))

輸出畫面

Hermite 差除表
i z(i) f(i)  f(i-1,i) f(i-2,i-1,i)...

1 0.0000 1.0000
2 0.0000 1.0000  1.0000
3 1.0000 2.7183  1.7183
3 1.0000 2.7183  1.7183 0.7183
4 1.0000 2.7183  2.7183 1.0000 0.2817
5 2.0000 7.3891  4.6708 1.9525 0.4762 0.0973
6 2.0000 7.3891  7.3891 2.7183 0.7658 0.1448 0.0238

Xa=1.5 Hn(xa)=4.481125---Hermite方法
真實解 f(x)=exp(x)
Xa=1.5 f(xa)=4.4816890703380645---真實值
Hn(x)與fn(x)的誤差0.000564070338064937

Julia語言例題3-1 求sin(x)的微分

Julia語言例題3-1 求sin(x)的微分

向前差 f'(x) = ( f(x+h) - f(x)) / h
向後差 f'(x) = ( f(x) - f(x-h)) / h

中央差 f'(x) = 0.5 * ( f(x+h) - f(x-h)) / h

取 h = 0.01

源自於
http://edisonx.pixnet.net/blog/post/90838030-%5Bc%E8%AA%9E%E8%A8%80%E6%95%B8%E5%80%BC%E5%88%86%E6%9E%90%5D-%E4%B8%80%E9%9A%8E%E5%BE%AE%E5%88%86%E6%B3%95

1. 若微分用定義方式做，效果有限，若欲求之 eps 過小 (如精度達 1E-9, 1E-12 )，一些情況上特別容易發生除零之窘境。

2. 數值分析之叢書不知道對微分法說明完不完整，即使跳過不講大概也不會讓人太過意外。但較意外的是，非線性方程式求解中的「牛頓法」，微分函式不是用數值分析方式算出來的，而是事先用筆算算出來的。

從定義做

一般微積分課本裡對於微分之定義大概如此

f ' (x) = [ f(x+h) - f(x) ] / h , h -> 0

當然可微包含的三個條件就不再重述了。

我們假設 f(x) = x * x ，要求 f'(2) 之值為何，

若取 h = 0.01，按上述的定義計算 f'(x) = [ f(2.01) - f(2) ] / 0.01 = 4.01；

再取 h = 0.001，按上述的定義計算 f'(x) = [f(2.001) - f(2) ] / 0.001 = 4.001;

實際上和 f'(2) 之誤差，與 h 成正比。

前差近似、後差近似、央央近似

然而上述所使用的模型，是屬較差的微分模型，所使用的模型為「二點前差」微分法。

在圖形上取得兩點，分別為 x，得到 f(x)，與 x 往前看一個 h 之 f(x+h)，

而微分所使用的方法絕不只這項，還有「後差除」、「中央差除」等方式，

所使用公式約略如下。

意到誤差項，那很重要，若對於輸出之微分誤差精度有所要求，

它可拿來做為概估之 h 推算。根據上述之三種不同方法，這次我們以 f(x) = sin(x) ，取 h = 0.01，

using Printf

#=====================================================
Trigonometric and hyperbolic functions
All the standard trigonometric and hyperbolic functions are also defined:
sin cos tan cot sec csc
sinh cosh tanh coth sech csch
asin acos atan acot asec acsc
asinh acosh atanh acoth asech acsch
sinc cosc atan2
=====================================================#

function fx(x::Float64) #// 欲微分函數
return sin(x)
end

function FordDiff(x::Float64 , h::Float64, fx::typeof(func)) #//前差微分
a=fx(x+h)
b=fx(x)
return ((a-b) / h)
end

function BackDiff(x::Float64 , h::Float64, fx::typeof(func)) #//後差微分
return ((fx(x) - fx(x-h) ) / h )
end

function MidDiff(x::Float64 , h::Float64, fx::typeof(func)) #//中差微分
return (0.5 * ( fx(x+h) - fx(x-h) ) / h)
end

#=====================================================#
x=1.0
h=0.01
answer = cos(x) # // 答案
cal = FordDiff(x, h, func)
delta = cal - answer
s=@sprintf("FordDiff : %0.5lf, delta = %0.5lf\n",cal,delta)
println(s)

cal = BackDiff(x, h, func)
delta = cal - answer
s=@sprintf("BackDiff : %0.5lf, delta = %0.5lf\n",cal,delta)
println(s)

cal = MidDiff(x, h, func)
delta = cal - answer
s=@sprintf("MidDiff : %0.5lf, delta = %0.5lf\n",cal,delta)
println(s)

輸出畫面

FordDiff : 0.53609, delta = -0.00422

BackDiff : 0.54450, delta = 0.00420

MidDiff : 0.54029, delta = -0.00001

Julia語言範例1-13 已知不明分行物體, 經光測儀器之量測,印出Hermite 差除表

Julia語言範例1-13 已知不明分行物體, 經光測儀器之量測,印出Hermite 差除表

範例1-13 已知不明分行物體, 經光測儀器之量測其結果如下
=========================================
時間位置飛行速度
x f(x) f ' (x)
0 50.00 50.00
2 216.67 94.44
4 410.00 98.00
=========================================
印出Hermite 差除表
當xa=1.5時求 Hn(xa)之值

程式

using Printf

#======================================================

i z(i) f(i) f(i-1,i) f(i-2,i-1,i)...

0 0.0 50.0000

1 0.0 50.0000 50.0000

2 2.0 216.6700 83.3350 16.6675

3 2.0 216.6700 94.4400 5.5525 -5.5575

4 4.0 410.0000 96.6650 1.1125 -1.1100 1.1119

5 4.0 410.0000 98.0000 0.6675 -0.2225 0.2219 -0.2225

function [H]=hermite(X,x,f,fd);

m=length(x);

for i=1:m

z(2*i-1)=x(i);

Q(2*i-1,1)=f(i);

z(2*i)=x(i);

Q(2*i,1)=f(i);

Q(2*i,2)=fd(i);

if i~=1

Q(2*i-1,2)=(Q(2*i-1,1)-Q(2*i-2,1))/(z(2*i-1)-z(2*i-2));

end;

for i=3:2*m

for j=3:i

Q(i,j)=(Q(i,j-1)-Q(i-1,j-1))/(z(i)-z(i-j+1));

end

p=1;

H=Q(1,1);

for i=2:2*m

p=p.*(X-z(i-1));

H=H+p.*Q(i,i);

end;

======================================================#

println("Hermite 差除表")

x=[0.0 , 2.0 ,4.0]

f=[50.00 , 216.67 , 410.00]

ff=[50.00 , 94.44 , 98.00]

n=length(x)

z=[0.0 for i=1:2*n]

q= [[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ],

[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ],

[ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ]]

println("i\tz(i)\tf(i)\t\tf(i-1,i)\tf(i-2,i-1,i)...\n");

for i=1:n

z[2*i-1]=x[i]

z[2*i]=x[i]

q[2*i-1][1]=f[i]

q[2*i][1]=f[i]

q[2*i][2]=ff[i]

if (i!=1)

q[2*i-1][2]= (q[2*i-1][1]- q[2*i-2][1]) / (z[2*i-1]-z[2*i-2])

end

#println(i,z,q)

end

for i=1:n

if(i==1)

s=@sprintf("%d\t%0.4f\t%0.4f\n",i,z[i],q[i][1])

print(s)

else (i==2)

s=@sprintf("%d\t%0.4f\t%0.4f\t\t%0.4f\n",i,z[i],q[i][1],q[i][2])

print(s)

end

for i=3:2*n

s=@sprintf("%d\t%0.4f\t%0.4f\t%0.4f",i,z[i],q[i][1],q[i][2])

print(s)

for j=3:i

q[i][j]=(q[i][j-1]-q[i-1][j-1])/(z[i]-z[i-j+1])

s=@sprintf("\t%0.4f",q[i][j])

print(s)

end

println()

end

xa=1.0
p=1
H=q[1][1]

for i=2:2*n
p *= (xa-z[i-1])
H = H+p*q[i][i]

#println(q[i][i])
#println(p)
#println(H)
#println()

end
println()

println("Xa=",xa," Hn(xa)=" ,H, "---Hermite方法")
H1=(100*xa + (50.0/(xa+1)))
println("Xa=",xa," f(xa)=" ,H1 , "---真實值")
println("Hn(x)與fn(x)的誤差" ,abs(H1-H))

輸出畫面

Hermite 差除表
i z(i) f(i)  f(i-1,i) f(i-2,i-1,i)...

1 0.0000 50.0000
2 0.0000 50.0000  50.0000
3 2.0000 216.6700 83.3350
3 2.0000 216.6700 83.3350 16.6675
4 2.0000 216.6700 94.4400 5.5525 -5.5575
5 4.0000 410.0000 96.6650 1.1125 -1.1100 1.1119
6 4.0000 410.0000 98.0000 0.6675 -0.2225 0.2219 -0.2225

Xa=1.0 Hn(xa)=124.004375---Hermite方法
Xa=1.0 f(xa)=125.0---真實值
Hn(x)與fn(x)的誤差0.995625000000004

訂閱：文章 (Atom)

alex9ufo 聰明人求知心切

2019年2月28日星期四

Julia語言習題3-6請選擇適當的微分近似法完成下表：

Julia語言微分的近似值習題3-3 完成表中的 f'(x) 與 f''(x)

Julia語言例題3-3 不等距的函數f(x)的微分近似值先使用牛頓向前的內插多項式

2019年2月25日星期一

Julia語言利用3點微分方式求f'(x)

Julia語言例題3-1 利用向前差　向後差　中央差　求f(x)=sin(x)的值

Julia語言例題 EX2-8 定點回路法求非線性方程式 f(x)=exp(x) - 3x^2 = 0

Julia語言例題 EX2-9 定點回路法求非線性方程式 f(x)=1/5^x - x = 0

Julia語言範例2-5 非線性方程式 f(x)= (4x-7) /(x-1)的根

julia語言牛頓法(Newton - Raphson)以圖形概念主要是不停在取切線斜率。

Julia語言例題2-6 已知方程式 e^x + x^-2 + 2 cosx -6 利用正割法找出f(x)=0的根

非線性方程式求解 - 目錄

Julia語言範例2-3 使用Newton-Raphson理則解 x^3 + 4x^2 -10 =0 非線性方程式解

Julia語言範例2-4 使用Newton-Raphson理則解 cos(x) - x =0 非線性方程式解

Julia語言 EX2-2 f(x) = (0.0000115pow(x,2)+0.000014pow(x,1.5)-0.01962) 利用二分法求解---(a) 測試根的範圍

2019年2月24日星期日

Julia語言 EX2-2 f(x)=(0.0000115pow(x,2)+0.000014pow(x,1.5)-0.01962) 二分法求解

Julia語言 EX2-2 f(x) = (0.0000115pow(x,2)+0.000014pow(x,1.5)-0.01962) 利用二分法求解

Julia語言 EX2-1利用2分法求 x^3-4x+2 =0 在0.0到2.0之間的解

Julia語言習題1-7 已知下面4點座標求Hermite內插法的插除表並計算Hn(x) = [1.5 , 2.5 , 3.5 , 4.5]之值

Julia語言例題1-12 已知3點座標及導數 ,建立 Hermite 向後差除表

Julia語言例題3-1 求sin(x)的微分

Julia語言範例1-13 已知不明分行物體, 經光測儀器之量測,印出Hermite 差除表

113 學年度第 1 學期 RFID應用課程 Arduino程式

2019年2月28日 星期四

2019年2月25日 星期一

2019年2月24日 星期日

2019年2月28日星期四

2019年2月25日星期一

2019年2月24日星期日