フィルターのクリア

Matlab to solve finite difference for Ode

2 ビュー (過去 30 日間)
Ojonile Joseph
Ojonile Joseph 2021 年 7 月 7 日
回答済み: Viranch Patel 2021 年 7 月 9 日
% Finite difference method
% Approximate the solution of y"=(-2/x)y'+(2/x^2)y+ sin(lnx)/x^2
% for 1<=x<=2 with y(1)=1 and y(2)=2.
p = @(x) (-2/x);
q = @(x) (2/x^2);
r = @(x) (sin(log(x))/x^2);
aa = 1; bb = 2; alpha = 1; beta = 2; n=29;      % h = (bb-aa)/(n+1);   h=0.1
fprintf('   x           w   \n');
h = (bb-aa)/(n+1);
a = zeros(1,n+1);
b = zeros(1,n+1);
c = zeros(1,n+1);
d = zeros(1,n+1);
l = zeros(1,n+1);
u = zeros(1,n+1);
z = zeros(1,n+1);
w = zeros(1,n+1);
x = aa+h;
a(1) = 2+h^2*q(x);
b(1) = -1+0.5*h*p(x);
d(1) = -h^2*r(x)+(1+0.5*h*p(x))*alpha;
m = n-1;
for i = 2 : m
x = aa+i*h;
a(i) = 2+h^2*q(x);
b(i) = -1+0.5*h*p(x);
c(i) = -1-0.5*h*p(x);
d(i) = -h^2*r(x);
end
x = bb-h;
a(n) = 2+h^2*q(x);
c(n) = -1-0.5*h*p(x);
d(n) = -h^2*r(x)+(1-0.5*h*p(x))*beta;
l(1) = a(1);
u(1) = b(1)/a(1);
z(1) = d(1)/l(1);
for i = 2 : m
l(i) = a(i)-c(i)*u(i-1);
u(i) = b(i)/l(i);
z(i) = (d(i)-c(i)*z(i-1))/l(i);
end
l(n) = a(n)-c(n)*u(n-1);
z(n) = (d(n)-c(n)*z(n-1))/l(n);
w(n) = z(n);
for j = 1 : m
i = n-j;
w(i) = z(i)-u(i)*w(i+1);
end
i = 0;
fprintf('%5.4f    %11.8f\n', aa, alpha);
for i = 1 : n
x = aa+i*h;
fprintf('%5.4f    %11.8f\n', x, w(i));
end
i = n+1;
fprintf('%5.4f    %11.8f\n', bb, beta);am trying to run d code but it keeps writing Error: Invalid text character. The text ' ' contains an unsupported non-ASCII whitespace character. Please where am I doing wrong?

回答 (1 件)

Viranch Patel
Viranch Patel 2021 年 7 月 9 日
I tried running the code on my machine and it's working fine. You can paste this code to your machine and see if its working. I am attaching the output too. You might be using some wrong character for any particular character which is not supported.
% Finite difference method
% Approximate the solution of y"=(-2/x)y'+(2/x^2)y+ sin(lnx)/x^2
% for 1<=x<=2 with y(1)=1 and y(2)=2.
p = @(x) (-2/x);
q = @(x) (2/x^2);
r = @(x) (sin(log(x))/x^2);
aa = 1; bb = 2; alpha = 1; beta = 2; n=29; % h = (bb-aa)/(n+1); h=0.1
fprintf(' x w \n');
h = (bb-aa)/(n+1);
a = zeros(1,n+1);
b = zeros(1,n+1);
c = zeros(1,n+1);
d = zeros(1,n+1);
l = zeros(1,n+1);
u = zeros(1,n+1);
z = zeros(1,n+1);
w = zeros(1,n+1);
x = aa+h;
a(1) = 2+h^2*q(x);
b(1) = -1+0.5*h*p(x);
d(1) = -h^2*r(x)+(1+0.5*h*p(x))*alpha;
m = n-1;
for i = 2 : m
x = aa+i*h;
a(i) = 2+h^2*q(x);
b(i) = -1+0.5*h*p(x);
c(i) = -1-0.5*h*p(x);
d(i) = -h^2*r(x);
end
x = bb-h;
a(n) = 2+h^2*q(x);
c(n) = -1-0.5*h*p(x);
d(n) = -h^2*r(x)+(1-0.5*h*p(x))*beta;
l(1) = a(1);
u(1) = b(1)/a(1);
z(1) = d(1)/l(1);
for i = 2 : m
l(i) = a(i)-c(i)*u(i-1);
u(i) = b(i)/l(i);
z(i) = (d(i)-c(i)*z(i-1))/l(i);
end
l(n) = a(n)-c(n)*u(n-1);
z(n) = (d(n)-c(n)*z(n-1))/l(n);
w(n) = z(n);
for j = 1 : m
i = n-j;
w(i) = z(i)-u(i)*w(i+1);
end
i = 0;
fprintf('%5.4f %11.8f\n', aa, alpha);
for i = 1 : n
x = aa+i*h;
fprintf('%5.4f %11.8f\n', x, w(i));
end
i = n+1;
fprintf('%5.4f %11.8f\n', bb, beta);
Output:
x w
1.0000 1.00000000
1.0333 1.03067949
1.0667 1.06155254
1.1000 1.09262608
1.1333 1.12390423
1.1667 1.15538887
1.2000 1.18708018
1.2333 1.21897697
1.2667 1.25107701
1.3000 1.28337728
1.3333 1.31587416
1.3667 1.34856355
1.4000 1.38144105
1.4333 1.41450201
1.4667 1.44774163
1.5000 1.48115505
1.5333 1.51473732
1.5667 1.54848354
1.6000 1.58238883
1.6333 1.61644834
1.6667 1.65065735
1.7000 1.68501118
1.7333 1.71950528
1.7667 1.75413522
1.8000 1.78889666
1.8333 1.82378540
1.8667 1.85879736
1.9000 1.89392857
1.9333 1.92917520
1.9667 1.96453354
2.0000 2.00000000

カテゴリ

Help Center および File ExchangeProgramming についてさらに検索

製品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by