對於那些看別人部落格,自己卻從來不寫的人,我希望你們也多寫。
對於那些只會說別人:''你寫的這個別人已經寫過了''的人,請你躝開。
****************************************==
用牛頓迭代法(newton-raphson)求解非線性方程組
x1**2-10*x1+x2**2+8=0
x1*x2**2+x1-10*x2+8=0
圖1是隨著求解次數的x1與x2值
圖2是隨著求解次數的誤差的二範數
求解非線性方程組
e1: x1**2-10*x1+x2**2+8=0
e2: x1*x2**2+x1-10*x2+8=0
寫出雅可比矩陣
e1對x1求導 e1對x2求導
e2對x1求導 e2對x2求導
即2*x1-10 2*x2
x2**2+1 2*x1*x2-10
牛頓迭代法 xn+1 = xn - f(xn)/f'(xn)
f1 = x1**2-10*x1+x2**2+8
f2 = x1*x2**2+x1-10*x2+8
df11 = 2*x1-10
df12 = 2*x2
df21 = x2**2+1
df22 = 2*x1*x2-10
"""import numpy as np
def calc(x):
"""傳入乙個向量x 然後求解得到f與df
"""x1=x[0,0] #x的0行0列 如果x[0]這樣引用的話 則x1與x2是乙個陣列而不是乙個數字
x2=x[1,0]
f1=x1**2-10*x1+x2**2+8
f2=x1*x2**2+x1-10*x2+8
df11=2*x1-10
df12=2*x2
df21=x2**2+1
df22=2*x1*x2-10
f=np.array([[f1,f2]]).t
df=np.array([[df11,df12],
[df21,df22]])
return f,df
#初始x1為0 x2為0
x0=np.array([[0,0]]).t
f0,df0=calc(x0)
#初始化
p_xn=x0
p_f=f0
p_df=df0
#定義乙個列表儲存迴圈過程中得到的所有x1 x2
x1 =
x2 =
error =
for i in range(100):
y = np.dot(np.linalg.inv(p_df),p_f)
p_xn1 = p_xn - y
p_xn = p_xn1
p_f,p_df = calc(p_xn)
#如果y的2範數小於0.000001則跳出迴圈
y_norm2 = np.linalg.norm(y,2)
if y_norm2<0.000001:
break
# 畫圖
import matplotlib.pyplot as plt
plt.figure(1)
plt.plot(x1,'r--',label='x1')
plt.plot(x2,'b-',label='x2')
plt.xlabel('solution times')
plt.legend()
plt.figure(2)
plt.plot(error,label='error')
plt.xlabel('solution times')
plt.legend()
c語言求解非線性方程組 C語言求解線性方程組
在之前的文章c語言實現矩陣求秩和化約化階梯形中,我們已經實現了求矩陣的秩與約化階梯形,在此基礎上,我們就可以來求解線性方程組了.一般線性方程組當且僅當 2.當 首先以增廣矩陣的形式輸入線性方程組,利用 matrix.h 標頭檔案中的求秩函式分別計算增廣矩陣和係數矩陣的秩,然後判斷是否有解 如果方程有...
MATLAB 線性方程組求解
clc,clear all close all 高斯消去法 a 2 3 4 3 5 2 4 3 30 線性方程組的係數矩陣 b 6 5 32 線性方程組的右端列向量 m,n size a 測量係數矩陣的維數 if m n fprint 線性方程組的係數矩陣非方陣 break end fprintf ...
MATLAB線性方程組求解
對於一般的,有唯一解的線性方程組,我們可以轉換成矩陣的形式 a x bax b ax b 則可以用矩陣運算求解x,即x a b 求解齊次線性方程組基礎解系的函式是null z null a 表示返回矩陣a的基礎解系組成的矩陣。z還滿足ztz i z null a,r 得出的z不滿足ztz i,但得出...