顯示具有 Julua程式語言 標籤的文章。顯示所有文章

2019年3月2日星期六

Julia語言例題4-6 雙重定積分 a=0.0 , b=1.0 , c=1.0 ,d=2.0 , f(x,y)= x^2 y 取n=10 求積分值

Julia語言例題4-6 雙重定積分 a=0.0 , b=1.0 , c=1.0 ,d=2.0 , f(x,y)= x^2 y 取n=10 求積分值

#==================================================
/* ex4-6-C.jl based on Trapezoidal Rule to
* compute the double integral.
*/
===================================================#
using Printf

function F(x::Float64, y::Float64) #// f(x,y)= x^2 * y
return (x*x*y)
end

function C(x::Float64) #// c(x)= 0.0
return (0.0)
end

function D(x::Float64) #// d(x)= 1.0
return (1.0)
end

function gy(n::Int64)
sum=0.0;
for i=0:n
for j=1:n-1
sum=sum+F(x[i+1],y[i+1][j+1])
#println(y[i][j])
end
g1[i+1]=(0.5*hy[i+1])*(F(x[i+1],y[i+1][1])+F(x[i+1],y[i+1][n+1])+2*sum)
println(i+1,"----",g1[i+1])
sum=0.0;
end
return g1
end

s=@sprintf("梯形積分計算雙重積分")
println(s)

x= [0.0 for i=1:20 ]
g1=[0.0 for i=1:20 ]
g2=[0.0 for i=1:20 ]
hy=[0.0 for i=1:20 ]

f(i) = [0.0 for i=1:20]
y= f.([0.0 for i=1:20])

sum=0.0
n=10
a=1.0
b=2.0

hx=(b-a)/n
ts=0.0
for i=0:n
x[i+1]=a+i*hx
hy[i+1]=(D(x[i+1])-C(x[i+1]))/n
for j=0:n
y[i+1][j+1]=C(x[i+1])+j*hy[i+1]
#println(y)
end
end

g2=gy(n)
sum1=0.0
println("\n\n")
for i=1:n-1
sum1=sum1+g2[i+1];
#println(g2[i+1])
end

ts= (hx/2) * (g1[1] + g1[n+1] + 2*sum1)
s=@sprintf("梯形積分計算雙重積分結果 T%d=%0.6lf\n",n,ts)
println(s)
s=@sprintf("實際值=%0.6lf\n",(7/6))
println(s)
tn=abs( 7/6 - ts )
s=@sprintf("誤差值=%0.6lf\n",tn )
println(s)

輸出畫面

梯形積分計算雙重積分
1----0.5000000000000001
2----0.6050000000000002
3----0.7200000000000001
4----0.8450000000000002
5----0.98
6----1.125
7----1.2800000000000005
8----1.4450000000000005
9----1.6200000000000003
10----1.8050000000000002
11----2.0000000000000004



梯形積分計算雙重積分結果 T10=1.167500

實際值=1.166667

誤差值=0.000833

Julia語言例題4-6 雙重定積分 a=0.0 , b=1.0 , c=1.0 ,d=2.0 , f(x,y)= x^2 y 取n=2 求積分值

Julia語言例題4-6 雙重定積分

b d

\ell =\iint _{T}f(x,y)\,dx\,dy

a c
a=0.0 , b=1.0 , c=1.0 ,d=2.0 , f(x,y)= x^2 y 取n=2 求積分值

#==================================================
/* ex4-6.jl based on Trapezoidal Rule to
* compute the double integral.
*/
===================================================#
using Printf

function F(x::Float64, y::Float64) #// f(x,y)= 4xy
return (x*x*y)
end

x=[0.0 for i=1:20]
g1=[0.0 for i=1:20 ]
y=[0.0 for i=1:20 ]

sum=0.0
n=2
a=1.0
b=2.0
hx=(b-a)/n
ax=a
for i=1:n+1
x[i]=ax
ax=ax+hx
end

c=0.0
d=1.0
hy=(d-c)/n
ay=c
for i=1:n+1
y[i]=ay
ay+=hy
end

for i=1:n+1
s=@sprintf("i=%2d , x[%1d]=%0.3f , y[%1d]=%0.3f ",i,i,x[i],i,y[i])
println(s)
end

println("\n梯形理則 (trapezoidal Rule) ")
Gx0=0.5*hy*(F(x[1],y[1]) + 2*F(x[1],y[2]) +F(x[1],y[3]))
s=@sprintf("Gx0=%0.5f",Gx0)
println(s)

Gx1=0.5*hy*(F(x[2],y[1])+2* F(x[2],y[2])+ F(x[2],y[3]) )
s=@sprintf("Gx1=%0.5f",Gx1)
println(s)

Gx2=0.5*hy*(F(x[3],y[1])+2* F(x[3],y[2])+ F(x[3],y[3]) )
s=@sprintf("Gx2=%0.5f",Gx2)
println(s)

Tn=0.5*hx*(Gx0+2*Gx1+Gx2)
s=@sprintf("0.5*hy*(Gx0+2*Gx1+Gx2)= %0.5f",Tn)
println(s)

#===============================================#
println("\n辛普森理則（Simpson's Rule)")
Gx0=(1/3)*hy*( F(x[1],y[1])+4* F(x[1],y[2])+ F(x[1],y[3]) )
s=@sprintf("Gx0=%0.5f",Gx0)
println(s)

Gx1=(1/3)*hy*( F(x[2],y[1])+4* F(x[2],y[2])+ F(x[2],y[3]) )
s=@sprintf("Gx1=%0.5f",Gx1)
println(s)

Gx2=(1/3)*hy*( F(x[3],y[1])+4* F(x[3],y[2])+ F(x[3],y[3]) )
s=@sprintf("Gx2=%0.5f",Gx2)
println(s)
Tn=(1/3)*hy*(Gx0+4*Gx1+Gx2)
s=@sprintf("(1/3)*hy*(Gx0+4*Gx1+Gx2)= %0.5f", Tn)
println(s)
#===============================================#

輸出畫面

i= 1 , x[1]=1.000 , y[1]=0.000 
i= 2 , x[2]=1.500 , y[2]=0.500 
i= 3 , x[3]=2.000 , y[3]=1.000 

梯形理則 (trapezoidal Rule) 
Gx0=0.50000
Gx1=1.12500
Gx2=2.00000
0.5*hy*(Gx0+2*Gx1+Gx2)= 1.18750

辛普森理則（Simpson's Rule)
Gx0=0.50000
Gx1=1.12500
Gx2=2.00000
(1/3)*hy*(Gx0+4*Gx1+Gx2)= 1.16667

2019年3月1日星期五

Julia語言例題4-4使用求 f(x)= (1/x) a=1.0 , b=2.0 求n=?? 辛普森積分法 , 梯形積分法需要幾次

Julia語言例題4-4使用
求 f(x)= (1/x) a=1.0 , b=2.0 求n=??

辛普森積分法需要 n 次 n=8
梯形積分法需要 n 次 n=41

x2
為了方便，符號簡化: ∫ dx = ∫ , h = x2-x0
x0
Thm:

Given [x0,x2] , x1 = (x0+x2)/2, assume f in C^4[x0,x2]

∫f = (h/3) * ( f(x0) + 4f(x1) + f(x2) ) - (h^5/90)*f^(4)(ξ)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~
這項是標準的辛普森法這項是誤差

#========================================================
/* ex4-4.jl based on Simpson's Rule to compute
* definite integral with domain [a,b] and
* n even-grid. n must be even.
*/
========================================================#
using Printf

function F4(x::Float64) #// 四次微分後函數
return (24/x^5)
end

function F1(x::Float64) #// 欲微分函數
return (1/x)
end

function F2(x::Float64) #// 二次微分後函數
return (2/x^3)
end

a=1.0
b=2.0
TOL=0.0001
xa=F4(a)
xb=F4(b)
if xa>xb
k=xa
else
k=xb
end

n1= k* (b-a)^5
n2= n1 / (180*TOL)
n3= sqrt(n2)
n4=sqrt(n3)
n=Int32(round(n4+0.5))
if mod(n,2)!=0
n=n+1
end
println("---------------------------------")
println("辛普森積分法需要",n,"次")

#驗證
println("---------------------------------")
sn=0.0
a=1.0
b=2.0
m=n/2
h=(b-a)/n
sum1=0.0
sum2=0.0
for i=1:2*m-1
x=a+i*h;
if(i%2==0)
sum2=sum2+F1(x)
sn= (h/3.0)*(F1(a)+F1(b)+2.0*sum2+4.0*sum1)

#s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F1(x))
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.6f err=%0.6f ",i,x,F1(abs(x)), abs(sn-(log(2)-log(1))))
println(s)
else
sum1=sum1+F1(x)
sn= (h/3.0)*(F1(a)+F1(b)+2.0*sum2+4.0*sum1)

#s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F1(x))
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.6f err=%0.6f ",i,x,F1(abs(x)), abs(sn-(log(2)-log(1))))
println(s)
end
end
print("辛普森積分法 ")
s=@sprintf("S%d=%lf 誤差=%0.6f\n",n,sn ,abs(sn-(log(2)-log(1))))
println(s)

#梯形法

a=1.0
b=2.0
TOL=0.0001
xa=F2(a)
xb=F2(b)
if xa>xb
k=xa
else
k=xb
end

n1= k* (b-a)^5
n2= n1 / (12*TOL)
n4= sqrt(n2)
n=Int32(round(n4+0.5))
println("---------------------------------")
println("梯形積分法需要",n,"次")

println("---------------------------------")
sn=0.0
a=1.0
b=2.0
h=(b-a)/n
x=a
result=0.0
for i=1:n-1
x=x+h
result=result+F1(abs(x))
tn=(h/2.0)*(F1(abs(a))+F1(abs(b))+2.0*result)
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.6f err=%0.6f ",i,x,F1(abs(x)), abs(tn-(log(2)-log(1))))
println(s)

end

#tn=(h/2.0)*(F1(abs(a))+F1(abs(b))+2.0*result)

print("梯形積分法 ")
s=@sprintf("S%d=%lf 誤差=%0.6f\n",n,tn ,abs(tn-(log(2)-log(1))))
println(s)

輸出畫面

---------------------------------
辛普森積分法需要8次
---------------------------------
i= 1 , x=1.13 ---- f(x)=0.888889 err=0.482499 
i= 2 , x=1.25 ---- f(x)=0.800000 err=0.415832 
i= 3 , x=1.38 ---- f(x)=0.727273 err=0.294620 
i= 4 , x=1.50 ---- f(x)=0.666667 err=0.239065 
i= 5 , x=1.63 ---- f(x)=0.615385 err=0.136501 
i= 6 , x=1.75 ---- f(x)=0.571429 err=0.088882 
i= 7 , x=1.88 ---- f(x)=0.533333 err=0.000007 
辛普森積分法  S8=0.693155 誤差=0.000007

---------------------------------
梯形積分法需要41次
---------------------------------
i= 1 , x=1.02 ---- f(x)=0.976190 err=0.651045 
i= 2 , x=1.05 ---- f(x)=0.953488 err=0.627789 
i= 3 , x=1.07 ---- f(x)=0.931818 err=0.605062 
i= 4 , x=1.10 ---- f(x)=0.911111 err=0.582840 
i= 5 , x=1.12 ---- f(x)=0.891304 err=0.561101 
i= 6 , x=1.15 ---- f(x)=0.872340 err=0.539824 
i= 7 , x=1.17 ---- f(x)=0.854167 err=0.518991 
i= 8 , x=1.20 ---- f(x)=0.836735 err=0.498582 
i= 9 , x=1.22 ---- f(x)=0.820000 err=0.478582 
i=10 , x=1.24 ---- f(x)=0.803922 err=0.458975 
i=11 , x=1.27 ---- f(x)=0.788462 err=0.439744 
i=12 , x=1.29 ---- f(x)=0.773585 err=0.420876 
i=13 , x=1.32 ---- f(x)=0.759259 err=0.402357 
i=14 , x=1.34 ---- f(x)=0.745455 err=0.384176 
i=15 , x=1.37 ---- f(x)=0.732143 err=0.366318 
i=16 , x=1.39 ---- f(x)=0.719298 err=0.348775 
i=17 , x=1.41 ---- f(x)=0.706897 err=0.331533 
i=18 , x=1.44 ---- f(x)=0.694915 err=0.314584 
i=19 , x=1.46 ---- f(x)=0.683333 err=0.297917 
i=20 , x=1.49 ---- f(x)=0.672131 err=0.281524 
i=21 , x=1.51 ---- f(x)=0.661290 err=0.265395 
i=22 , x=1.54 ---- f(x)=0.650794 err=0.249522 
i=23 , x=1.56 ---- f(x)=0.640625 err=0.233897 
i=24 , x=1.59 ---- f(x)=0.630769 err=0.218512 
i=25 , x=1.61 ---- f(x)=0.621212 err=0.203361 
i=26 , x=1.63 ---- f(x)=0.611940 err=0.188435 
i=27 , x=1.66 ---- f(x)=0.602941 err=0.173729 
i=28 , x=1.68 ---- f(x)=0.594203 err=0.159237 
i=29 , x=1.71 ---- f(x)=0.585714 err=0.144951 
i=30 , x=1.73 ---- f(x)=0.577465 err=0.130867 
i=31 , x=1.76 ---- f(x)=0.569444 err=0.116978 
i=32 , x=1.78 ---- f(x)=0.561644 err=0.103279 
i=33 , x=1.80 ---- f(x)=0.554054 err=0.089765 
i=34 , x=1.83 ---- f(x)=0.546667 err=0.076432 
i=35 , x=1.85 ---- f(x)=0.539474 err=0.063274 
i=36 , x=1.88 ---- f(x)=0.532468 err=0.050287 
i=37 , x=1.90 ---- f(x)=0.525641 err=0.037467 
i=38 , x=1.93 ---- f(x)=0.518987 err=0.024809 
i=39 , x=1.95 ---- f(x)=0.512500 err=0.012309 
i=40 , x=1.98 ---- f(x)=0.506173 err=0.000037 
梯形積分法  S41=0.693184 誤差=0.000037

Julia語言例題4-4利用辛普森理則（Simpson's Rule)與梯形法計算定積分

Julia語言例題4-4利用辛普森理則（Simpson's Rule) 與梯形法
# (a) 計算 exp( 1/x) 在[1, 2]的定積分
# (b) 計算 ( 1/x) 在[1, 2]的定積分

#========================================================
/* ex4-4.jl based on Simpson's Rule to compute
* definite integral with domain [a,b] and
* n even-grid. n must be even.
*/
========================================================#
using Printf

function F1(x::Float64) #// 欲微分函數
return (exp(1/x))
end

function F2(x::Float64) #// 欲微分函數
return (1/x)
end

a=1.0
b=2.0
n=4
m=n/2
h=(b-a)/n
sum1=0.0
sum2=0.0
for i=1:2*m-1
x=a+i*h;
if(i%2==0)
sum2=sum2+F1(x)
s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F1(x))
println(s)
else
sum1=sum1+F1(x)
s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F1(x))
println(s)
end
end
sn= (h/3.0)*(F1(a)+F1(b)+2.0*sum2+4.0*sum1)
print("辛普森積分法 ")
s=@sprintf("S%d=%lf\n",n,sn)
println(s)

println("---------------------------------")
n=10
a=1.0
b=2.0
n=4
h=(b-a)/n
x=a
result=0.0
for i=1:n-1
x=x+h
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.6f",i,x,F1(abs(x)))
println(s)
result=result+F1(abs(x))
end

tn=(h/2.0)*(F1(abs(a))+F1(abs(b))+2.0*result)

print("梯形積分法 ")
s=@sprintf("T%d=%0.7lf\n",n,tn)
println(s)

println("---------------------------------")
n=10
m=n/2
h=(b-a)/n
sum1=0.0
sum2=0.0
for i=1:2*m-1
x=a+i*h;
if(i%2==0)
sum2=sum2+F2(x)
s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F2(x))
println(s)
else
sum1=sum1+F2(x)
s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F2(x))
println(s)
end
end
sn= (h/3.0)*(F2(a)+F2(b)+2.0*sum2+4.0*sum1)
print("辛普森積分法 ")
s=@sprintf("S%d=%lf\n",n,sn)
println(s)

println("---------------------------------")
n=10
a=1.0
b=2.0
n=10
h=(b-a)/n
x=a
result=0.0
for i=1:n-1
x=x+h
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.6f",i,x,F2(abs(x)))
println(s)
result=result+F2(abs(x))
end

tn=(h/2.0)*(F2(abs(a))+F2(abs(b))+2.0*result)

print("梯形積分法 ")
s=@sprintf("T%d=%0.7lf\n",n,tn)
println(s)

輸出畫面

i= 1 ,x=1.250 --- F(x)=2.225541
i= 2 ,x=1.500 --- F(x)=1.947734
i= 3 ,x=1.750 --- F(x)=1.770795
辛普森積分法  S4=2.020651

---------------------------------
i= 1 , x=1.25 ---- f(x)=2.225541
i= 2 , x=1.50 ---- f(x)=1.947734
i= 3 , x=1.75 ---- f(x)=1.770795
梯形積分法  T4=2.0318929

---------------------------------
i= 1 ,x=1.100 --- F(x)=0.909091
i= 2 ,x=1.200 --- F(x)=0.833333
i= 3 ,x=1.300 --- F(x)=0.769231
i= 4 ,x=1.400 --- F(x)=0.714286
i= 5 ,x=1.500 --- F(x)=0.666667
i= 6 ,x=1.600 --- F(x)=0.625000
i= 7 ,x=1.700 --- F(x)=0.588235
i= 8 ,x=1.800 --- F(x)=0.555556
i= 9 ,x=1.900 --- F(x)=0.526316
辛普森積分法  S10=0.693150

---------------------------------
i= 1 , x=1.10 ---- f(x)=0.909091
i= 2 , x=1.20 ---- f(x)=0.833333
i= 3 , x=1.30 ---- f(x)=0.769231
i= 4 , x=1.40 ---- f(x)=0.714286
i= 5 , x=1.50 ---- f(x)=0.666667
i= 6 , x=1.60 ---- f(x)=0.625000
i= 7 , x=1.70 ---- f(x)=0.588235
i= 8 , x=1.80 ---- f(x)=0.555556
i= 9 , x=1.90 ---- f(x)=0.526316
梯形積分法  T10=0.6937714

Julia語言例題4-5 利用辛普森理則（Simpson's Rule) 計算 x=0 to 2 , 以曲線 y= (1+x^3) ^(1/3) 沿 z 軸旋轉一周的體積

Julia語言例題4-5 利用辛普森理則（Simpson's Rule)
#計算 x=0 to 2 , 以曲線 y= (1+x^3) ^(1/3) 沿 z 軸旋轉一周的體積

#========================================================
/* ex4-5.jl based on Simpson's Rule to compute
* definite integral with domain [a,b] and
* n even-grid. n must be even.
*/
========================================================#
using Printf

function F(x::Float64) #// 欲微分函數
return (( (1+x^3)^ (1/3) )^2)

end
a=0.0
b=2.0
n=10
m=n/2
h=(b-a)/n
sum1=0.0
sum2=0.0
for i=1:2*m-1
x=a+i*h;
if(i%2==0)
sum2=sum2+F(x)
s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F(x))
println(s)
else
sum1=sum1+F(x)
s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F(x))
println(s)
end
end
sn= pi *(h/3.0)*(F(a)+F(b)+2.0*sum2+4.0*sum1)
print("\n\n辛普森積分法 ")
s=@sprintf("S%d=%lf\n",n,sn)
println(s)

輸出畫面

i= 1 ,x=0.200 --- F(x)=1.005326
i= 2 ,x=0.400 --- F(x)=1.042224
i= 3 ,x=0.600 --- F(x)=1.139259
i= 4 ,x=0.800 --- F(x)=1.317350
i= 5 ,x=1.000 --- F(x)=1.587401
i= 6 ,x=1.200 --- F(x)=1.952374
i= 7 ,x=1.400 --- F(x)=2.411148
i= 8 ,x=1.600 --- F(x)=2.961326
i= 9 ,x=1.800 --- F(x)=3.600520


辛普森積分法  S10=12.325078

Julia語言例題EX4-2 請用梯形法求統計學常態分布的面積

f(x) = \frac{1}{\sqrt{2\pi}} \, \exp\left(-\frac{x^2}{2} \right)

。

a=-3.5 b=3.0 n=10 的積分面積

程式
#================================================
* ex4-2.jl based on Trapezoidol Rule is
* used for computing definite integral with
* domain [a,b] with n even-grids.
*
梯形法 T'= (h/2) [ f(a) + f(b) + 2 sigma i=1 to n-1 f(xi)]

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

using Printf

function f(x::Float64) #// 欲微分函數
return return (1.0/exp(x*x/2))
end

result=0.0
n=10
a=-3.5
b=-3.0

h=(b-a)/n
x=a
for i=1:n-1
x=x+h
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.4f",i,x,f(abs(x)))
println(s)
result=result+f(abs(x))
end

tn=(1.0/sqrt(2* pi ))*(h/2.0)*(f(abs(a))+f(abs(b))+2.0*result)

println("\nf(x)=(1.0/exp(x*x/2)) a=0.0 ,b=1.0 n=10 ")
s=@sprintf("T%d=%0.7lf\n",n,tn)
println(s)

ET1= (b^2 -1)* exp(- (b^2)/2)
#println(ET1)
ET= (ET1 * (b-a)^3) / (12*n^2)
s=@sprintf("梯形法的誤差範圍= %0.7f",ET)
println(s)

輸出畫面

i= 1 , x=-3.45 ---- f(x)=0.0026
i= 2 , x=-3.40 ---- f(x)=0.0031
i= 3 , x=-3.35 ---- f(x)=0.0037
i= 4 , x=-3.30 ---- f(x)=0.0043
i= 5 , x=-3.25 ---- f(x)=0.0051
i= 6 , x=-3.20 ---- f(x)=0.0060
i= 7 , x=-3.15 ---- f(x)=0.0070
i= 8 , x=-3.10 ---- f(x)=0.0082
i= 9 , x=-3.05 ---- f(x)=0.0095

f(x)=(1.0/exp(x*x/2)) a=0.0 ,b=1.0 n=10 
T10=0.0011194

梯形法的誤差範圍= 0.0000093

Julia語言利用梯形積分法求f (x) = exp(x) a=0.0 , b=1.0 n=10的積分值

Julia語言例題4-1利用梯形積分法求f (x) = exp(x) a=0.0 , b=1.0 n=10的積分值

#================================================
* ex4-1.jl based on Trapezoidol Rule is
* used for computing definite integral with
* domain [a,b] with n even-grids.
*
梯形法 T'= (h/2) [ f(a) + f(b) + 2 sigma i=1 to n-1 f(xi)]

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

using Printf

function f(x::Float64) #// 欲微分函數
return exp(x)
end

result=0.0
n=10
a=0.0
b=1.0

h=(b-a)/n
x=a
for i=1:n-1
x=x+h;
s=@sprintf("x=%0.1f ---- f(x)=%0.4f",x,f(abs(x)))
println(s)
result=result+f(abs(x))
end

tn=(h/2.0)*( f(abs(a)) + f(abs(b)) + 2.0*result )

println("\nf(x)=exp(x) a=0.0 ,b=1.0 n=10 ")
s=@sprintf("T%d=%10.6lf\n",n,tn)
println(s)

s=@sprintf("梯形積分與實際值的誤差值=%0.5f",abs(tn-(f(1.0)-f(0.0))))
println(s)

ET= (f(b)* (b-a)^3) / (12*n^2)
s=@sprintf("梯形法的誤差範圍= %0.5f",ET)
println(s)

輸出畫面

x=0.1 ---- f(x)=1.1052
x=0.2 ---- f(x)=1.2214
x=0.3 ---- f(x)=1.3499
x=0.4 ---- f(x)=1.4918
x=0.5 ---- f(x)=1.6487
x=0.6 ---- f(x)=1.8221
x=0.7 ---- f(x)=2.0138
x=0.8 ---- f(x)=2.2255
x=0.9 ---- f(x)=2.4596

f(x)=exp(x) a=0.0 ,b=1.0 n=10 
T10=  1.719713

梯形積分與實際值的誤差值=0.00143

梯形法的誤差範圍= 0.00227

Julia語言例題4-3 梯形法求標準常態分配的面積

f(x) = {1 \over \sigma\sqrt{2\pi} }\,e^{- {{(x-\mu )^2 \over 2\sigma^2}}}

f(x) = \frac{1}{\sqrt{2\pi}} \, \exp\left(-\frac{x^2}{2} \right)

a=-5 , b=5 , n=100

#================================================
* ex4-3.jl based on Trapezoidol Rule is
* used for computing definite integral with
* domain [a,b] with n even-grids.
*
梯形法 T'= (h/2) [ f(a) + f(b) + 2 sigma i=1 to n-1 f(xi)]

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

using Printf

function f(x::Float64) #// 欲微分函數
return (1.0/exp(x*x/2))
end

result=0.0
n=100
a=-5.0
b=5.0

h=(b-a)/n
x=a
for i=1:n
x=x+h;
result=result+f(abs(x))
end

tn=(1.0/sqrt(2* pi ))*(h/2.0)*(f(abs(a))+f(abs(b))+2.0*result)

s=@sprintf("T%d=%10.6lf\n",n,tn)
println(s)

輸出畫面

T100=  1.000000

Julia語言習題3-4請使用中央差近似法計算一次微分 f '(1.005) , f ' (1.015) 與二次微分 f '' (1.01) 之值

Julia語言習題3-4 觀察下列數據
x = 1.00 , 1.01 , 1.02
f(x)=1.27 , 1.32 , 1.38
請使用中央差近似法計算
一次微分 f '(1.005) , f ' (1.015) 與二次微分 f '' (1.01) 之值

程式
#中央近似差方法
# f' i = (f(i+1) - f(i-1)) / 2h
# f'' i = ( fi+1 - 2 fi + f i-1) / h^2
#============================================
x f(x) f'(x) f''(x)
1.00 1.27
1.005
1.01 1.32
1.015
1.02 1.38
============================================#

using Printf

x=[ 1.00 , 1.01 , 1.02 ]
f= [ [1.27 ,0.0 , 0.0 ],
[1.32 ,0.0 , 0.0 ],
[1.38 ,0.0 , 0.0 ] ]

xa=[1.005 , 1.015 ,1.01]
h=abs(x[2]-x[1])
i=1
s=@sprintf("f(%0.3f)的微分.......",xa[i])
println(s)
fd1=(f[i+1][1]- f[i][1]) / (h)
s=@sprintf("一次微分值= %0.4f ",fd1)
f[i][2]=fd1
println(s)

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

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

println("\n")
i=2
s=@sprintf("f(%0.3f)的二次微分.......",xa[i+1])
println(s)

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

輸出畫面

f(1.005)的微分.......
一次微分值= 5.0000 


f(1.015)的微分.......
一次微分值= 6.0000 


f(1.010)的二次微分.......
二次微分值= 100.0000

Julia語言習題3-3 等距的函數f(x)的微分近似值

Julia語言習題3-3 等距的函數f(x)的微分近似值

先使用牛頓向前的內插多項式
x= f(x)= f ' (x)= f ''(x)=
=============================================
0.15 0.1761 2.8775 -17 前向
0.17 0.2304 2.5675 -14.75 中央公式相同
0.19 0.2788 2.295 -12.5 中央公式相同
0.21 0.3222 2.06 -10.25 後向
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

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

訂閱：文章 (Atom)

alex9ufo 聰明人求知心切

2019年3月2日星期六

Julia語言例題4-6 雙重定積分 a=0.0 , b=1.0 , c=1.0 ,d=2.0 , f(x,y)= x^2 y 取n=10 求積分值

Julia語言例題4-6 雙重定積分 a=0.0 , b=1.0 , c=1.0 ,d=2.0 , f(x,y)= x^2 y 取n=2 求積分值

2019年3月1日星期五

Julia語言例題4-4使用求 f(x)= (1/x) a=1.0 , b=2.0 求n=?? 辛普森積分法 , 梯形積分法需要幾次

Julia語言例題4-4利用辛普森理則（Simpson's Rule)與梯形法計算定積分

Julia語言例題4-5 利用辛普森理則（Simpson's Rule) 計算 x=0 to 2 , 以曲線 y= (1+x^3) ^(1/3) 沿 z 軸旋轉一周的體積

Julia語言例題EX4-2 請用梯形法求統計學常態分布的面積

Julia語言利用梯形積分法求f (x) = exp(x) a=0.0 , b=1.0 n=10的積分值

Julia語言例題4-3 梯形法求標準常態分配的面積

Julia語言習題3-4請使用中央差近似法計算一次微分 f '(1.005) , f ' (1.015) 與二次微分 f '' (1.01) 之值

Julia語言習題3-3 等距的函數f(x)的微分近似值

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的根

MQTT 協定與 Modbus 通訊的遠端監控與控制系統架構

2019年3月2日 星期六

2019年3月1日 星期五

2019年2月28日 星期四

2019年2月25日 星期一

2019年3月2日星期六

2019年3月1日星期五

2019年2月28日星期四

2019年2月25日星期一