求线性方程组的系数

373965956

对于以上的线性方程，我知道ε11、ε22、ε33、γ12、γ13、γ23、σij的值，想求E、G和泊松比，在matlab中怎么编程啊，谢谢了


注：由于是初次在matlab板块中发帖，有不符合的地方，请版主帮忙修改

lengyunfeng

如果你的stress和strain的数值是理论对应值的话，可以用下面的方法计算。若为实测值的话，该法提供的解可能会因为你的测量精度而有一定的测量误差和计算误差，仅供参考。

1. 
   
2. %提问者UID：327380
   
3. %提问时间：2010-9-10
   
4. %主题名：求线性方程组的系数
   
5. %提问领域：符号方程
   
6. function [E,U,G]=wuxingjie(stress,strain)
   
7. %E:弹模，MPa
   
8. %U:泊松比
   
9. %G：剪切模量
   
10. %stress：应力向量（1*6或者6*1），MPa，前三个为各面上的正应力，后三个为剪应力
    
11. %strain：应变向量（1*6或者6*1），前三个为各面上的正应变，后三个为剪应变
    
12. %运行环境：winxp+matlab 7.0
    
13. %将stress与strain向量转变化列向量，以便以矩阵操作
    
14. if size(stress,2)==6
    
15.     stress=stress';
    
16. end
    
17. if size(strain,2)==6
    
18.     strain=strain';
    
19. end
    
20. syms e u g%定义三个未知量为字符变量
    
21. a=(eye(3)-ones(3,3)+diag([1/u,1/u,1/u]))*u/e;%生成系数阵左上角子阵
    
22. b=eye(3)/g;%生成系数阵右下角子阵
    
23. c=[a,zeros(3,3);zeros(3,3),b];%生成系数阵
    
24. d=c*stress-strain;%构建方程对应的函数
    
25. [e1,u1]=solve(d(1),d(2));%解e和u
    
26. g1=solve(d(4));%解g
    
27. %将三个字符精确解转变为双精度解
    
28. E=double(e1);
    
29. U=double(u1);
    
30. G=double(g1);
    

复制代码

373965956

谢谢斑竹的回答，我用这个代码计算了，可是27，28，29行出错，可能是字符型转换到数值型时出现错误

lengyunfeng

3# 373965956
错误提示帖出来嘛，我这没出错。而且我觉得如果真的你的机子上运行出错的话，应该也不会提示27、28、29三行出错吧？matlab运行到一句出错它就会停下来不往下走的

scott198510

2# lengyunfeng

请问版主这样一个数据拟合怎么求解呢？

比如这样一个公式：
H=A0*(F^n)*exp(-Q/(RT));  %  A0, n 为待求的参数
Q=275000;R=8.31;  %两个参数为定值

现在有几组数据：即T取确定值时，对应如下：

T=1073 时,   F=50*（1e+6) ,  H=1e-3
T=1053  时，F=68*(1e+6)    ,H=1e-3 ;
                    
T=1023  时，F=72*(1e+6)   , H=1e-3 ;
                     
T=973  时， F=120*(1e+6)    ,  H=1e-3 ;


                    


现在需要采用线性回归法求出待定的参数，
就是说求出的参数代入原始方程，各组数据尽可能吻合
哪位指点一下：

lengyunfeng

5# scott198510

你可以把四种情况组合成6个方程组，然后分别求解A0和n，再求这6种情况的平均值。下面的代码供你参考。

1. [a12,n12]=solve('1e-3=x*(50*1e+6)^y*exp(-275000/(8.31*1073))',...
   
2.     '1e-3=x*(68*1e+6)^y*exp(-275000/(8.31*1053))');
   
3. [a13,n13]=solve('1e-3=x*(50*1e+6)^y*exp(-275000/(8.31*1073))',...
   
4.     '1e-3=x*(72*1e+6)^y*exp(-275000/(8.31*1023))');
   
5. [a14,n14]=solve('1e-3=x*(50*1e+6)^y*exp(-275000/(8.31*1073))',...
   
6.     '1e-3=x*(120*1e+6)^y*exp(-275000/(8.31*973))');
   
7. [a23,n23]=solve('1e-3=x*(68*1e+6)^y*exp(-275000/(8.31*1053))',...
   
8.     '1e-3=x*(72*1e+6)^y*exp(-275000/(8.31*1023))');
   
9. [a24,n24]=solve('1e-3=x*(68*1e+6)^y*exp(-275000/(8.31*1053))',...
   
10.     '1e-3=x*(120*1e+6)^y*exp(-275000/(8.31*973))');
    
11. [a34,n34]=solve('1e-3=x*(72*1e+6)^y*exp(-275000/(8.31*1023))',...
    
12.     '1e-3=x*(120*1e+6)^y*exp(-275000/(8.31*973))');
    
13. A0=(a12+a13+a14+a23+a24+a34)/6;
    
14. N=(n12+n13+n14+n23+n24+n34)/6;
    
15. clear a** n**
    

复制代码

求出来的A0=.88917143988272878205794635502172e-5;N=5.5978087442400244774059171249202

lin2009

拟合之前数据要处理一下。
１、Ｈ当成常数；
２、120*(1e+6)不是很规范的写法；
３、曲线拟合用１stopt求解很方便，请转到MathCad / 1stopt 版块。

lin2009

6# lengyunfeng

不是平均值，是最小二乘解。

qibbxxt

我觉得可以用非线性拟合

1. function yhat = myTestFun(beta,x)
   
2. A0 = beta(1);
   
3. n = beta(2);
   
4. T=x(:,1);
   
5. F=x(:,2);
   
6. Q=275000;
   
7. R=8.31;
   
8. yhat = A0*(F.^n).*exp(-Q./(R*T));

复制代码

1. clear;clc;close all
   
2. x=[1073,50*10^6;1052,68*10^6;1023,72*10^6;973,120*10^6];
   
3. y=10^(-3)+zeros(size(x,1),1);
   
4. % beta=[0.89;5.59];
   
5. beta=[1;1]
   
6. beta = nlinfit(x,y,@myTestFun,beta)

复制代码

根据物理意义去调整初值，应该可以得到较好的结果

lin2009

参数        最佳估算
----------        -------------
a0        4.42504642359754E-20
n        3.84754755953478

====== 结果输出 ======

目标t        计算t
1073        1073.00000000542
1053        1052.99999998543
1023        1023.00000000695
973        973.000000002794

====== 计算结束 ======
迭代数: 60
计算用时(时:分:秒:毫秒): 00:00:18:531
计算中止原因: 达到收敛判定标准
优化算法: 通用全局优化法(UGO1)
目标函数值(最小): 3.06612790329988E-7
均方差(RMSE): 8.62827076231963E-9
残差和(RSS): 2.97788225391599E-16
相关系数(R): 1
决定系数(DC): 1

相应的1stopt程序如下：
Constant Q = 275000;
Constant R = 8.31;                   //% 两个参数为定值
Constant H = 0.001;                 // H直接作为常数，可降低拟合方程的维数。
Parameter A0,n;                      // 待拟合的参数
Variable  F,T;                          // 拟合变量
Function  H = A0*(F^n)*exp(-Q/(R*T)) ;//拟合方程
Data;
50000000, 1073
68000000, 1053
72000000, 1023
120000000, 973

scott198510

10# lin2009
我想问问楼主大虾，你这种算法是在1stopt里面做出来的，在MATLAB里面可以直接用类似的程序来非线性拟合吗？我不太会用1stopt,因为实际要处理的数据也就是像前面的方程有很多个，数据用的越多才会使得拟合的系数结果越准确

qibbxxt

11# scott198510
10#是将问题转化为优化问题求解了
请看9#的程序，就是用matlab作非线性优化的程序

373965956

4# lengyunfeng
对不起，前几天太忙，现在才来回版主的帖子，我把代码贴出来

1. function [E,U,G]=wuxingjie(stress,strain)
   
2. %E:弹模，MPa
   
3. %U:泊松比
   
4. %G：剪切模量
   
5. %stress：应力向量（1*6或者6*1），MPa，前三个为各面上的正应力，后三个为剪应力
   
6. %strain：应变向量（1*6或者6*1），前三个为各面上的正应变，后三个为剪应变
   
7. %运行环境：winxp+matlab 7.0
   
8. %将stress与strain向量转变化列向量，以便以矩阵操作
   
9. stress=[20.0298E-06     20.0298E-06    -145.785E-06    -14.9050E-21    -12.1764E-06    -12.1764E-06];
   
10. strain=[-4.52486E+06    -4.52486E+06    -45.6361E+06    -1.84773E-09    -1.50947E+06    -1.50947E+06];
    
11. if size(stress,2)==6
    
12.     stress=stress';
    
13. end
    
14. if size(strain,2)==6
    
15.     strain=strain';
    
16. end
    
17. syms e u g%定义三个未知量为字符变量
    
18. a=(eye(3)-ones(3,3)+diag([1/u,1/u,1/u]))*u/e;%生成系数阵左上角子阵
    
19. b=eye(3)/g;%生成系数阵右下角子阵
    
20. c=[a,zeros(3,3);zeros(3,3),b];%生成系数阵
    
21. d=c*stress-strain;%构建方程对应的函数
    
22. [e1,u1]=solve(d(1),d(2));%解e和u
    
23. g1=solve(d(4));%解g
    
24. %将三个字符精确解转变为双精度解
    
25. E=double(e1);
    
26. U=double(u1);
    
27. G=double(g1);

复制代码

以下是出错提示

1. ??? Error using ==> mupadmex
   
2. Error in MuPAD command: DOUBLE cannot convert the input expression into a double array.
   
3. 
   
4. If the input expression contains a symbolic variable, use the VPA function instead.
   
5. 
   
6. Error in ==> sym.sym>sym.double at 936
   
7.             Xstr = mupadmex('symobj::double', S.s, 0);
   
8. 
   
9. Error in ==> strain_stress at 25
   
10. E=double(e1);
    

复制代码

我的matlab版本是2010.
如果我把后三行改成：fprintf('%d %d %d','e1 u1 g1')就可以得到正确答案。

scott198510

12# qibbxxt

请问版主，为什么运行你提供的代码时出现错误，而且出现不了所求的系数的值

qibbxxt

14# scott198510
我的代码不会出现错误，只是你运行的方法不对
第一段是函数文件，你保存为myTestFun.m
然后在command window里面输入第二段，或者建立脚本运行
结果的好与不好和初值有关系

lengyunfeng

13# 373965956
被你雷到了。。。你是不是把我给你的程序保存成 strain_stress这个文件名了？看你的提示第9行里说的是Error in ==> strain_stress at 25，估计错就错在这了；scott198510运行qibbxxt的程序出错的原因也在这，肯定是你把他第一段程序保存成乱七八糟的名字了

scott198510

用斑竹的方法求出的数据拟合的还比较好，因为用前面的数据求出来的与实际偏差太大，我把数据改成这样的：
就是说在温度确定的情况下:T=1053
H=1e-4,   F=24*（10^6）;
H=3e-4,   F=58*（10^6）;
H=1e-3,   F=70*（10^6）;
H=3e-3,   F=107*（10^6）;

参照斑竹的方法，把原来的表达式做一下变形，Q=160000值改小了一点，做出来的程序是这样的，可是为什么求出来的数据变成了虚数呢？即使把数据取绝对值，反带入方程里面求出的数据也感觉拟合的不好，远没有之前的数据那样对每组数据都比较接近：

function That = myTestFun(beta,x)
A0 = beta(1);
n = beta(2);
yb=x(:,1);
F=x(:,2);
Q=160000;
R=8.31;
That = Q/R./log(A0.*(F.^n)./yb);


x=[1e-4,24*10^6;3e-4,58*10^6;1e-3,70*10^6;3e-3,107*10^6];
y=1053+zeros(size(x,1),1);
% beta=[0.89;5.59];
beta=[1;1];
beta = nlinfit(x,y,@myTestFun,beta)；
A0 = beta(1)
n = beta(2)



求出来的数据是这样的：
A0 =
  8.8590e-008 +1.3118e-007i

n =
   1.4836 - 0.0546i

373965956

16# lengyunfeng
首先感谢qibbxxt和lengyunfeng两位版主的回答，呵呵，看了你们的解答，我也汗啊，我这是第一次用matlab，呵呵。看来还要多多学习啊

scott198510

9# qibbxxt
请问斑竹为什么用你的方法拟合的效果是这样的：

A0 =
  9.2511e-014

n =
    3.0418

这与上面楼主用优化解法求出
A0        4.42504642359754E-20
n        3.84754755953478
差的这麽远呢？而实际上用优化解法得到的数据要更加准确

scott198510

10# lin2009

用楼主的结果反代入方程发现拟合的非常好，而且与实际的参考数据
也相等吻合，想问问楼主可不可以帮忙再看看下面我17#回复的帖子里面的问题怎么求法，小弟对mathcad不懂，只会使用MATLAB，但是MATLAB用非线性方法解出的数据差的太远了