Matlab to solve finite difference for Ode

Question

Ojonile Joseph 2021-7-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/873958-matlab-to-solve-finite-difference-for-ode

回答： 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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Viranch Patel 2021-7-9

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/873958-matlab-to-solve-finite-difference-for-ode#answer_743478

在 MATLAB Online 中打开

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   
0000     1.00000000
0333     1.03067949
0667     1.06155254
1000     1.09262608
1333     1.12390423
1667     1.15538887
2000     1.18708018
2333     1.21897697
2667     1.25107701
3000     1.28337728
3333     1.31587416
3667     1.34856355
4000     1.38144105
4333     1.41450201
4667     1.44774163
5000     1.48115505
5333     1.51473732
5667     1.54848354
6000     1.58238883
6333     1.61644834
6667     1.65065735
7000     1.68501118
7333     1.71950528
7667     1.75413522
8000     1.78889666
8333     1.82378540
8667     1.85879736
9000     1.89392857
9333     1.92917520
9667     1.96453354
0000     2.00000000