解非線性方程式

[shyfang (水深火熱)]

%in Matlab

% 想要解出x1,x2的方程式, a1,a2,a3 和y1,y2是data
%        a1*x1 + a2*x1*x2 = y1
%        a3*x2 + a2*x1*x2 = y2
%
% i.e 使f1=0 f2=0 的(x1,x2)
%        f1 = a1*x1 + a2*x1*x2 - y1
%        f2 = a3*x2 + a2*x1*x2 - y2
%
%
function f = fun(x,a,y)
        f=zeros(2,1);
        f(1)=a(1)*x(1)+a(2)*x(1)*x(2)-y(1);
        f(2)=a(3)*x(2)+a(2)*x(1)*x(2)-y(2);

end
%通常是把function那一段另外存一個*.m檔，放在同一個目錄下

y=[ 2.1,0.5];
a=[200199000000,58546439234,1];
%解出0根
h=fsolve(@(x) fun(x,a,y),[0;0]);


# in R
# 雖然用Matlab解很方便，但是matlab 2007非常佔資源
# 在R裡面要用optim (找最小值的function)
# function改成找||f||的。使得norm最小就行了
# 想解的只有x，是x的函式
fun=function(x){
        f=rep(0,length(x));
        f[1]=a[1]*x[1]+a[2]*x[1]*x[2]-y[1]
        f[2]=a[3]*x[2]+a[2]*x[1]*x[2]-y[2]
        normf=sum(f^2)           # f 的 norm^2
        return(normf)
}

y=c( 2.1,0.5);
a=c(200199000000,58546439234,1);

# 把 y&a 的值 assign 好了之後就可以使用
opt<-optim(rep(0,nunknown),fun)

# 因為function不能丟其它的變數，要小心資料的污染。

--
※ 發信站: 批踢踢實業坊(ptt.cc) 
◆ From: 140.114.36.142
推 catherinejoy:有個 stanford 作的 cvx 可在matlab 上使用 也很不賴