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

10.8.1 打靶法及其MATLAB程序

10.8.1  打靶法及其MATLAB程序

用线性打靶法和常用的四阶龙格-库塔公式(10.75)求解线性边值问题(10.95)数值解的MATLAB主程序

   

function [k,X,Y,wucha,P]=xxdb(dydx1,dydx2,a,b,alpha,beta,h)

n=fix((b-a)/h); X=zeros(n+1,1); CT1=[alpha,0];

Y=zeros(n+1,length(CT1)); Y1=zeros(n+1,length(CT1));

Y2=zeros(n+1,length(CT1));

X=a:h:b;

Y1(1,:)= CT1;

CT2=[0,1];Y2(1,:)= CT2;

for k=1:n

k1=feval(dydx1,X(k),Y1(k,:))

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

k2=feval(dydx1,x2,y2);

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

k4=feval(dydx1, X(k)+h,Y1(k,:)'+k3*h);

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

end

u=Y1(:,1)

for k=1:n

k1=feval(dydx2,X(k),Y2(k,:))

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

k2=feval(dydx2,x2,y2);

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

k4=feval(dydx2, X(k)+h,Y2(k,:)'+k3*h);

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

end

v=Y2(:,1)

Y=u+(beta-u(n+1))*v/v(n+1)

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'];

plot(X,Y(:,1),'ro',X,Y1(:,1),'g*',X,Y2(:,1),'mp')

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

legend('是边值问题的数值解y(x)的曲线','是初值问题1的数值解u(x)的曲线', '是初值问题2的数值解v(x)的曲线')

title('用线性打靶法求线性边值问题的数值解的图形')

 

例10.8.1  用线性打靶法求线性边值问题

, 的数值解,h=0.5,并与精确解比较,画出它们的图形.

解  在MATLAB工作窗口输入程序

>> a=0;b=6;h=0.5; alpha=1.25;beta=-1.25; c1=2;c2=2;c3=0.8;

[k,X,Y,wucha,P]=xxdb(@dydx1,@dydx2,a,b,alpha,beta,h),

hold on

y=1.2088219648217098859769672950885*X+1.25-2.05*X.^2+8/5*X.*atan(X)-2/5*log(1+X.^2)+2/5*log(1+X.^2).*X.^2;

plot(X,y,'b-'), hold off

legend('是边值问题的数值解Y(x)的图形','是初值问题1的数值解u(x)的图形', '是初值问题2的数值解v(x)的图形','是边值问题的精确解y(x)的曲线')

title('用线性打靶法求线性边值问题的数值解和精确解的图形')

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

for k=1:n+1

wuchay(k)=norm(y(k)-Y(k)); k=k+1;

end

wuchay;

k=1:n+1;wuchay=wuchay(1:k,:);

P1=[k',X',y',Y, wuchay' ,wucha']

运行后求得此线性边值问题在自变量X 处的数值解的值Y, 偏差,精确解

的值y 及其与精确解的绝对误差,它们的图形.

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