what can i further do in my this code of linear shooting
2 ビュー (過去 30 日間)
古いコメントを表示
function[]=LS(x,y) clc % Dy is the representation of y' %D2y is the representation od y'' D2y=inline('-(2/x)*Dy + (2/x^2)*y + sin(ln(x))/x^2','x','y','Dy') p=inline('-(2/x)','x'); p1=p; q=inline('(2/x^2)','x'); q1=q; r=inline('sin(log(x))/x^2','x'); r1=r; a=input('enter the lower limit of x: '); b=input('enter the upper limit of x: '); N=input('enter the integer N: '); y0=input('enter the lower boundary value of y: '); y1=input('enter the upper boundary value of y: '); fprintf('The system of linear Differential Equations the become\n') yA=inline('-(2/x)*Dy + (2/x^2)*y + sin(ln(x))/x^2','x','y','Dy') fprintf('with the conditions y(%1g)=%1g and Dy(%1g)=0\n',a,y0,a) f1=yA; yB=inline('-(2/x)*Dy + (2/x^2)*y','x','y','Dy') fprintf('with the conditions y(%1g)=0 and Dy(%1g)=1\n',a,a) f2=yB; h=(b-a)/N; u1=y0; u2=0; v1=0; v2=1; w1=y0; fprintf('\n x u1 v1 w1\n'); fprintf('%10g %10g %10g %10g\n',a,u1,v1,w1); for i=0:N-1; x=a+i*h; K11=h*(u2); K12=h*(((p1(x)*u2))+(q1(x)*u1)+r1(x)); K21=h*(u2+(1/2)*K12); K22=h*((p1(x+h/2)*(u2+(1/2)*K12)+q(x+h/2)*(u1+(1/2)*K11)+r(x+h/2))); K31=h*(u2+(1/2)*K22); K32=h*((p1(x+h/2)*(u2+(1/2)*K22)+q(x+h/2)*(u1+(1/2)*K21)+r(x+h/2))); K41=h*(u2+K32); K42=h*((p(x+h)*(u2+K32)+(q(x+h)*(u1+K31))+r(x+h))); u1=u1+(1/6)*(K11+2*K21+2*K31+K41); u2=u2+(1/6)*(K12+2*K22+2*K32+K42);
k11=h*(v2);
k12=h*(((p1(x)*v2))+(q1(x)*v1));
k21=h*(v2+(1/2)*k12);
k22=h*((p1(x+h/2)*(v2+(1/2)*k12))+(q(x+h/2)*(v1+(1/2)*k11)));
k31=h*(v2+(1/2)*k22);
k32=h*((p1(x+h/2)*(v2+(1/2)*k22))+(q(x+h/2)*(v1+(1/2)*k21)));
k41=h*(v2+k32);
k42=h*((p(x+h)*(v2+k32))+(q(x+h)*(v1+k31)));
v1=v1+(1/6)*(k11+2*k21+2*k31+k41);
v2=v2+(1/6)*(k12+2*k22+2*k32+k42);
fprintf('%10g %10g %10g\n',x+h,u1,v1);
end
4 件のコメント
回答 (1 件)
Jan
2011 年 7 月 4 日
Do not use "k24" if you want to access an array - "k(2, 4)" is smarter, faster and the musltiplication by "h" can be performed for all elements simultaneoulsly with "h*k".
0 件のコメント
参考
カテゴリ
Help Center および File Exchange で Waveform Generation についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!