matlab数值分析程序--高等数学,数值代数的matlab实现-文字版, matlab电子书, 和matlab 有关的电子书:

10.7.1 一阶微分方程组的数值解及其MATLAB程序

10.7.1  一阶微分方程组的数值解及其MATLAB程序

 

例10.7.2  求微分方程组在区间H=[0.1,60]上满足条件:时,的特解.

解  调用dzdx2.m求解,在MATLAB工作窗口输入程序

>> H=[0.1,60];

z0=[1;1;1];

[x,z]=ode15s('dzdx2',H,z0);

plot(x,z(:,1),'g-',x,z(:,2),'b*',x,z(:,3),'mp')

xlabel('轴\it x'); ylabel('轴\it y')

legend('是方程解z1的曲线','是方程解z2的曲线', '是方程解z3的曲线')

运行后求得解函数 ,的图形.

 

2、根据求微分方程组的数值解的公式编写MATLAB程序计算

用常用的四阶龙格-库塔公式(10.75)求解常微分方程组的初值问题(10.67)的数值解的MATLAB主程序

function [k,X,Y,wucha,P]=RK4z(dydx,a,b,CT,h)

n=fix((b-a)/h);

X=zeros(n+1,1);

Y=zeros(n+1,length(CT));

X=a:h:b;

Y(1,:)= CT';

for k=1:n

k1=feval(dydx,X(k),Y(k,:))

x2=X(k)+h/2;y2=Y(k,:)'+k1*h/2;

k2=feval(dydx,x2,y2);

k3=feval(dydx,x2,Y(k,:)'+k2*h/2);

k4=feval(dydx, X(k)+h,Y(k,:)'+k3*h);

Y(k+1,:)=Y(k,:)+h*(k1'+2*k2'+2*k3'+k4')/6;k=k+1;

end

for k=2:n+1

wucha(k)=norm(Y(k)-Y(k-1)); k=k+1;

end

X=X(1:n+1);Y=Y(1:n+1,:);k=1:n+1;wucha=wucha(1:k,:);

P=[k',X',Y,wucha'];

 

例10.7.3  用常用的四阶龙格-库塔公式(10.75)求解例10.7.2的常微分方程组的初值问题的数值解,h=0.25,画出数值解的图形,并与例10.7.2比较.

解 ⑴ 建立名为dydx.m的M函数文件

function  dY=dydx(X,Y)

dY(1)=Y(2);

dY(2)=Y(3);

dY(3)=Y(3)*X^(-1)-3*X^(-2)*Y(2)+2*X^(-3)*Y(1)+9*X^3*sin(X);

dY=[dY(1);dY(2);dY(3)];

⑵ 调用dydx.m求解,在MATLAB工作窗口输入程序

>> CT=[1;1;1];h=0.25;

[k,X,Y,wucha,P]=RK4z(@dydx,0.1,60,CT,h),

plot(X,Y(:,1),'g-',X,Y(:,2),'b*',X,Y(:,3),'mp')

xlabel('轴\it x'); ylabel('轴\it y')

legend('是方程解z1的曲线','是方程解z2的曲线', '是方程解z3的曲线')

运行后求得解函数 ,的数值解(略)及其图形.

欢迎转载,转载请注明来自一手册:http://yishouce.com/book/3/3100701.html
友情链接It题库(ittiku.com)| 版权归yishouce.com所有| 友链等可联系 admin#yishouce.com|粤ICP备16001685号-1