声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1109|回复: 6

随便问问

[复制链接]
发表于 2006-5-19 15:38 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
老师要求用动画表示(最好) 界面也可以 关于数值分析的<BR>好难 下面是代码 有高手帮忙看下知道不<BR>function [x,det,index] = Gauss(A,b)<BR>% 求线性方程组的列主元 Gauss 消去法,其中,<BR>% A 为方程组的系数矩阵;<BR>% b 为方程组的右端项;<BR>% x 为方程组的解;<BR>% det 为系数矩阵A的行列式的值;<BR>% index 为指标变量,index=0表示计算失败,index=1表示计算成功。<BR>[n,m] = size(A);nb  = length(b);<BR>% 当方程组行与列的维数不相等时,停止计算,并输出出错信息。<BR>if n~=m <BR>    error('The rows and columns of matrix A must be equal!');<BR>    return;<BR>end<BR>%当方程组与右端项的维数不匹配时,停止计算,并输出出错信息<BR>if m~=nb<BR>    error('The columns of A must be equal the length of b!');<BR>    return;<BR>end<BR>%开始计算,先赋初值<BR>index = 1; det = 1;x=zeros(n,1);<BR>for k=1:n-1<BR>    % 选主元<BR>    a_max = 0;<BR>    for i = k:n<BR>        if abs(A(i,k))&gt;a_max<BR>            a_max=abs(A(i,k));r=1;<BR>        end<BR>    end<BR>    if a_max&lt;1e-10<BR>        index = 0;return;<BR>    end <BR>    % 交换两行<BR>    if r&gt;k<BR>        for j=k:n<BR>            z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;<BR>        end<BR>        z=b(k);b(k)=b(r);b(r)=z;det=-det;<BR>    end<BR>    %消元过程<BR>    for i=k+1:n<BR>        m=A(i,k)/A(k,k);<BR>        for j = k+1:n<BR>            A(i,j)=A(i,j)-m*A(k,j);<BR>        end<BR>        b(i)=b(i)-m*b(k);<BR>    end<BR>    det=det*A(k,k);<BR>end<BR>det=det*A(n,n);<BR>% 回代过程。<BR>if abs(A(n,n))&lt;1e-10<BR>    index=0;return;<BR>end<BR>for k=n:-1:1<BR>    for j=k+1:n<BR>        b(k)=b(k)-A(k,j)*x(j);<BR>    end<BR>    x(k)=b(k)/A(k,k);<BR>end
回复
分享到:

使用道具 举报

发表于 2006-5-19 18:42 | 显示全部楼层

回复:(sys1983)随便问问

什么意思?你想问什么问题?
 楼主| 发表于 2006-5-20 09:09 | 显示全部楼层
<P>教授你好 对上面的 线性方程组求解的过程 加点代码<BR> 我们老师要求用动画表示(最好), 如或者做个界面也可以 输入矩阵A, b 后面得到结果x的解 </P>
发表于 2006-5-20 14:04 | 显示全部楼层

回复:(sys1983)随便问问

是不是上述函数永动画的形式画图?<BR><BR>这个前面有一个摆线的例子,你搜索一下,把里边的函数换成你的函数就基本差不多了,可能个别参数需要调整一下
 楼主| 发表于 2006-5-20 14:19 | 显示全部楼层
教授再问一下 也是数值分析的<BR>引用了MATLAB的函数<STRONG>ode23<BR></STRONG>这个列向量为yprime。同样,y1和y2合并写成列向量y。所得函数M文件是:<BR>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">function yprime=vdpol(t , y);</FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">%  VDPOL(t , y) returns derivatives of the Van der Pol equation:</FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">%  </FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">%  x </FONT>‘‘<FONT face="Times New Roman">-mu *(1-x ^2)*x </FONT>‘<FONT face="Times New Roman">+x=0 (</FONT>‘<FONT face="Times New Roman"> = d/dx , </FONT>‘‘<FONT face="Times New Roman"> = d^2/dx^2)</FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">%  </FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">%  let y(1)=x and y(2)=x</FONT>‘</P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">%  </FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">%  then y(1)</FONT> ‘<FONT face="Times New Roman"> = y(2)</FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">%      y(2)</FONT> ‘<FONT face="Times New Roman"> = MU*(1-y(1)^2)*y(2)-y(1)</FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><p><FONT face="Times New Roman"> </FONT></p></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">global MU %   choose 0&lt;MU&lt;10 in Command workspace</FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><p><FONT face="Times New Roman"> </FONT></p></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">yprime=[y(2)  MU*(1-y(1)^2)*y(2)-y(1)];  %  output must be a column</FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><p><FONT face="Times New Roman"> </FONT></p></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt">计算结果如下:</P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><p><FONT face="Times New Roman"> </FONT></p></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">&gt;&gt;global MU %  define MU as a global variable in the Command Workspace</FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">&gt;&gt;MU=2;  %  set global parameter to desired value</FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">&gt;&gt;[t , y]=ode23(‘</FONT> <FONT face="Times New Roman">vdpol ’</FONT><FONT face="Times New Roman"> , 0 , 30 , [1 ; 0]);  %  to=0 , tf=30 , yo=[1 ; 0]</FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">&gt;&gt;y1=y( : , 1);  %  first column is y(1) versus time points in t</FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">&gt;&gt;y2=y( : , 2);  %  second column is y(2)</FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">&gt;&gt;plot(t , y1 , t , y2 , </FONT>‘ <FONT face="Times New Roman">-- </FONT>‘<FONT face="Times New Roman">)</FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">&gt;&gt;xlabel(</FONT>‘<FONT face="Times New Roman"> Time , Second </FONT>‘<FONT face="Times New Roman">) , ylabel(</FONT>‘<FONT face="Times New Roman"> Y(1) and Y(2) </FONT>‘<FONT face="Times New Roman">)</FONT></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">&gt;&gt;title(</FONT>‘<FONT face="Times New Roman"> Van der Pol Solution for mu=2 </FONT>‘<FONT face="Times New Roman">)<BR>前面是对的 按例应该可以出图的 (后面有问题)<BR></P>
<P 0cm 0cm 0pt; TEXT-INDENT: 24pt"><FONT face="Times New Roman">[t , y]=ode23(</FONT>‘ <FONT face="Times New Roman">vdpol </FONT>‘<FONT face="Times New Roman"> , 0 , 30 , [1 ; 0]);  %  to=0 , tf=30 , yo=[1 ; 0]<BR>我替换成了<BR>       <FONT face="Times New Roman">[t , y]=ode23(‘</FONT><FONT face="Times New Roman">vdpol</FONT><FONT face="Times New Roman"> ’, [0 , 30] , [1 , 0]);  <BR>但是后面还是有错    <BR>能帮我看下指导下吗?</FONT></FONT></P></FONT>
发表于 2006-5-20 14:27 | 显示全部楼层

回复:(sys1983)随便问问

<P 24pt? TEXT-INDENT: 0pt; 0cm><FONT face="Times New Roman">函数改为<BR><BR>function yprime=vdpol(t , y);<BR>global MU %   choose 0&lt;MU&lt;10 in Command workspace<BR>yprime=[y(2); MU*(1-y(1)^2)*y(2)-y(1)];  %  output must be a column<BR><BR>调用改为<BR>[t , y]=ode23('vdpol' , [0,30] , [1,0]);<BR><BR><BR></FONT></P>
 楼主| 发表于 2006-5-21 11:04 | 显示全部楼层
<P>好的 谢谢教授 </P>
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-29 17:56 , Processed in 0.073490 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表